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SYNCHRONOUS  DEMODULATION  OF 
CHOPPER  MODULATED  RADIANCE  SIGNALS 

Initiator:  Mr.  T.  Murdock 

Problem  No. :  3003 
Project  No. :  7670 

This  problem  entailed  the  development  of  a  sequence  of  programs  to  proc¬ 
ess  telemetered  radiance  signals.  The  main  processing  was  performed  in 
a  program  called  IRDSP.  The  following  contains  a  description  of  the  proc¬ 
essing  algorithms  developed  by  RDP  for  this  task. 

Program  IRDSP  (Infra- Red- Data  Signal  Processing)  was  written  to  de¬ 
modulate  the  chopper-modulated  radiance  signals  contained  in  the  tele¬ 
metered  data  received  from  the  IRBS  or  ZIP  flights.  The  multi-tape  telem¬ 
etry  data  base  is  processed  by  IRDSP  one  tape  at  a  time.  The  output  of 
the  program  consists  of  plots  and  printed  tables  of  the  radiance  estimates 
and  a  plot -data  file.  When  the  entire  input  data  base  has  been  processed, 
the  plot -data  files  are  sorted  by  key  and  a  plot-data  base  is  created.  Using 
this,  flight-plots  of  the  radiance  signals  are  produced  using  a  separate  pro¬ 
gram.  The  flight  plots  can  be  over  the  entire  flight  time,  or  over  a  time 
interval  falling  on  two  input  tapes. 

In  this  write-up  the  description  of  program  IRDSP  is  presented  in  four 
sections: 

1.  The  Signal  Demodulation  task. 

2.  Input  Data  Base. 

3.  Processing  Methodology  and  Program  Flow  Chart. 

4.  Demodulation  Algorithm. 
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1.  The  Signal  Demodulation  Task 

The  high  speed  data  link  in  both  the  IRBS  and  ZIP  instrumentation  transmit 
the  radiance  measurements  made  by  the  multi-sensor  Infra-Red  (IR)  de¬ 
tectors  in  a  chopper  modulated  form.  27  of  the  38  channels  in  IRBS  and  32 
of  the  44  channels  in  ZIP  are  chopper-modulated  outputs  of  Infra  Red  (IR) 
detectors.  Each  is  sensitive  to  a  unique  range  of  wavelengths.  RDP's  task 
is  to  demodulate  the  radiance  message  in  each  channel  and  deliver  plots 
of  time  versus  radiance  tables  to  the  client. 

A  typical  chopper -modulated  radiance  waveform  in  the  input  data  base 
has  the  form  as  shown  in  Figure  1.  The  radiance  message  is  embedded  in 
the  relatively  slowly  varying  envelope  of  the  chopping  waveform. 

Over  the  20-minute  flight  time,  the  detectors  are  exposed  to  both 
calibrated  and  unknown  target  sources  of  radiance.  Each  channel  output 
will  contain  the  numbered  features  shown  in  Figure  1.  These  are  briefly 
described  below. 

1.  Just  prior  to  take-off  the  detectors  are  calibrated  by  making  the  hard¬ 
ware  execute  an  auto-calibration  (Autocal)  sequence  in  which  known 
levels  of  black  body  radiation  are  shone  on  the  detectors. 

2.  Approximately  60  seconds  after  take-off,  the  autocal  sequence  is  re¬ 
peated  for  an  inflight  calibration. 

3.  After  100  seconds  into  the  flight,  the  detectors  scan  the  earth  limb  in 
1/2°  increments.  This  data  is  taken  for  the  next  200  secs. 

4.  During  the  remainder  of  the  flight,  the  detectors  measure  Zodiacal  light. 
To  design  a  suitable  demodulation  algorithm,  we  must  have  an  adequate 

mathematical  representation  of  the  received  chopped  waveforms. 

The  chopping  of  the  radiance  message  is  done  by  mechanically  opening 
and  closing  the  detector  aperture  in  step  with  a  high-Q  tuning  fork  of  fre¬ 
quency  fc  Hz. 

The  chopping  waveform  is  not  rectangular.  Nevertheless  it  is  periodic 
and  can  be  represented  by  a  Fourier  series  as: 

»  jnw  t 

<P(t)  *  ^2  an  ’e  °  *  (1) 

n=-<*> 
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Its  Fourier  Transform  is  a  train  of  impulses. 


CP 

4>(w)  =  2 ff  •  6(w-nwc)  . 


(2) 


n=-oo 


When  the  radiance  message,  r(t)  is  chopped  by  <p(t),  the  convolution  theo¬ 
rem  states  that  r(t)  •  <p(t)  has  the  spectrum 


F(r(t)  -<p(t))  =^R(w)*  4>(w) 


(where  F(  )  is  the  Fourier  Transform  operator  and  *  is  the  convolution 
operator) 


=  R(w)  * 


w  -nw  ) 
c 


•  R(w  -nwc) 


(3) 


i.e, ,  the  chopped  message  spectrum  is  the  weighted  sum  of  the  translates 
of  the  message  spectrum,  each  translate  being  located  at  a  chopper  har¬ 
monic. 

Before  sampling  and  PCM  encoding,  however,  the  message  spectrum 
at  D.C.  is  removed  by  a  band-pass  filter  which  also  acts  as  a  de-aliasing 
filter.  Now,  if  the  message  spectrum  R(w)  has  a  bandwidth  B«wc,  then, 


s(t)  ~  r(t)  •  <p (t ) I  <4> 

I  Band  Pass  Filtered 

has  the  approximate  spectra 

QO 

s(w)  =  ^  an  '  Gn  '  R<w-nwc)  ’  {5) 

n=-a> 

where  G  are  the  values  of  the  filter’s  frequency  response  at  the  n  har- 
n 

monic  location  and  is  assumed  to  be  constant  in  a  local  neighborhood, 

nf  ±  B  Hz, 
c 


4 


Letting  A^  =  a^  •  Gn>  the  received  chopped  radiance  signals  then  have 
spectra  of  the  form, 

CD 

S(w)  =  Y  An  *  R(w-nwc)  .  (6) 

n=-® 

Note  that  A  =  0  for  n  =  0  and  n  >  f  /2f  due  to  the  band-pass  filter, 
n  sc 

The  signal  processing  task  is  to  reconstruct  the  radiance  message  r(t) 
from  a  noisy  observation  of  s(t).  In  practice,  this  can  be  done  only  to  with¬ 
in  a  scale  factor  since  the  Fourier  coefficients,  an,  of  the  chopping  wave¬ 
form  as  manifested  in  s(t)  are  not  known.  The  proper  scale  factor  is  then 
determined  from  the  autocalibration  level  outputs. 

2.  Input  Data  Base 

The  ZIP  and  IRBS  data  bases  have  different  amounts  and  types  of  data,  but 
as  far  as  radiance  data  demodulation  is  concerned,  both  data  bases  have 
the  same  structure  and  are  thus  processed  in  the  same  manner.  The  IRBS 
rocket  is  expected  to  transmit  information  for  about  20  minutes.  The  on¬ 
board  electronic  hardware  samples  38  desired  output  quantities,  quantizes 
each  sample  to  14  bits,  multiplexes  and  block-encodes  the  binary  data  to 
facilitate  receiver  synchronization  and  transmits  the  output  bit  stream  at 
300  kilo-bits  per  sec  (kbps)  as  a  PCM  Biphase  (L)  signal.  At  the  receiver, 
the  received  bit  sequence  is  recorded  on  an  analog  tape.  In  the  ZIP  design 
44  output  quantities  are  sampled,  multiplexed,  and  telemetered  at  616  kbps. 
In  either  case,  the  analog  data  tape  is,  however,  not  directly  accessible  by 
a  computer.  In  fact,  to  make  the  data  compatible  with  program  IRDSP 
several  steps  are  involved.  These  steps  are  block  diagrammed  in  Figure  2. 
A  brief  description  of  the  steps  in  Figure  2  is  as  follows: 

1.  The  analog  tape  is  decoded  in  pieces  on  to  "n"  digital  tapes  by  the 
Honeywell  316  FM  decoder.  Each  tape  (designated  1. 1-I.n)  can  con¬ 
tain  approximately  3  minutes  of  IRBS  data  or  about  90  secs  of  ZIP 
data. 

2.  Duplication  on  to  work  tapes.  The  tapes  in  Step  1,  i.e.,  1. 1-I.n  are 
simply  duplicated  on  to  work  tapes,  II.  1-II.n.  All  further  steps  are 
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Figure  2.  Steps  in  the  Decoding  and  Restructuring  of  the  IRBS/ZIP  Analog  Data  Tape. 


performed  on  these  work  tapes.  The  reason  is  that  in  case  of  loss  of 
information  in  any  tape,  recovery  via  duplication  in  Step  2  has  much 
faster  turnaround  time  than  recovery  via  decoding  in  Step  1. 

3.  Unpacking/ Repacking.  The  data  in  tapes  II.  1-II.n  (the  outputs  from 
Step  2)  exists  in  blocks  of  packed  bits.  A  modified  version  of  HON316 
is  used  to  unpack  and  repack  specific  bits  into  integer  values.  These 
numbers  are,  in  fact,  the  output  of  the  A/D  converter  in  the  rocket 
hardware  which  originally  sampled  the  data  channels.  HON316  then 
writes  these  integer  sample  values  in  blocks  to  the  output  tape. 

Each  output  tape  (III.  1-III.n)  contains  one  file  comprised  of  binary 
records  of  the  form  shown  in  Figure  3. 


*1 

*2 

‘13s  IRBS 
t  :  ZIP 

dl“l> 

• 

• 

• 

• 

• 

...  *  .  ... 

d38(t1):  IRBS 

d  (t  ):  IRBS 

d44^11>:  ZIP 

Figure  3.  Structure  of  a  Record  In  the  Input  Tapes  for  Program 
IRDSP. 

Each  record  contains  the  two-integer  dimensions  of  a  data  array 
followed  by  the  array  itself.  Array  dimensions  are  38  X 13  for  IRBS 
and  44  xll  for  ZIP.  The  first  location  of  each  column  contains  sam¬ 
pling  time  values  and  the  remaining  locations  contain  channel  outputs 
at  this  time  instant. 

In  the  ideal  situation  of  no  time  errors  (row  #1  of  each  record), 
the  tapes  III.  1-III.n  are  compatible  with  program  IRDSP.  Practically, 
however,  errors  are  invariably  present  and  the  following  step  deals 
with  this  problem. 
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4.  Identification  and  classification  of  time  errors.  Program  ANALYZE 
is  used  to  read  the  tapes  III.  1— III.  n  (one  tape  per  run)  and  scan  through 
the  first  row  of  each  record  Identifying  and  classifying  time-data  er¬ 
rors.  It  also  produces  a  plot  of  time  versus  sample  index  in  which 
errors  can  be  visually  observed. 

Each  run's  printout  is  next  inspected  and  a  judgment  is  made  on 
whether  the  corresponding  tape  is  acceptable.  Tapes  that  are  deemed 
unacceptable  are  returned  to  be  re-digitized. 

The  above  four-step  procedure  for  preparing  the  raw,  rocket  data  for 
analysis  by  program  IRDSP,  appears  deceptively  straight  forward.  In  fact, 
it  is  a  time-consuming  process  due  to  the  turnaround  requirements  for  tape- 
jobs  on  the  AFGL  computer  system.  Experience  with  the  test  data  tape  has 
shown  that  4-5  weeks  are  needed  to  complete  Steps  1-4  above,  if  the  data 
exhibits  a  large  number  of  time  errors.  Further,  the  probability  exists 
that  the  time  errors  will  lead  to  misinterpretation  of  the  data. 

3.  Processing  Methodology  and  Program  Flowchart 
Earlier,  in  Section  2,  the  final  structure  of  the  received  rocket  data  was 
described.  The  data  exists  on  9-track  tapes.  Consider  the  IRBS  data  for 
a  moment.  Each  tape  has  approximately  3  minutes  of  data.  Since  the  sys¬ 
tem  sampling  rate  is  563.91  Hz,  one  tape  contains  over  100,000  samples 
of  data  for  each  of  the  38  channels.  For  processing  on  the  AFGL  computer 
which  has  a  300k  word  memory,  one  is  forced  to  process  the  data  in  smaller 
convenient  segments.  Currently,  a  segment  size  of  4000  is  being  used.  This 
is  approximately  7  seconds  of  data.  This  segmentation  does  not,  however, 
interfere  with  the  signal  processing.  Successive  input  segments  are  precisely 
overlapped  to  compensate  for  data  loss  due  to  end-effects  of  the  digital  fil¬ 
lers  and  consequently  the  output  blocks  are  contiguously  spaced  in  time. 

Of  course,  if  a  data  gap  is  detected,  the  segmenting  process  continues  from 
beyond  the  gap.  The  process  continues  till  the  data  corresponding  to  the 
specified  time  interval  has  been  analysed  or  until  an  end-of-file  mark  has 
been  reached. 

The  implementation  of  the  above  processing  methodology  is  in  the  form 

of  two  routines  -  FNDBLKS  and  GETBLK.  FNDBLKS  is  called  once  before 
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any  signal  processing  begins.  It  scans  through  the  input  data  tape  and  iden¬ 
tifies  segments  to  be  processed.  This  identification  involves  tagging  the 
location  in  the  file  of  the  beginning  time  and  end  time  for  each  segment  in 
terms  of  (1)  the  record  #  in  the  file  and  (2)  the  column  #  within  this  record. 

The  second  routine,  GETBLK,  uses  the  information  compiled  by 
FNDBLKS  to  read  a  data  block  from  tape,  convert  the  integer  A/D  output 
values  into  voltage  values  and  make  the  block  available  to  the  subroutine 
ENVELOP  for  signal  processing.  After  a  call  in  which  a  new  channel  is 
specified,  successive  calls  to  GETBLK  can  be  used  to  sequentially  access 
the  data  segments  for  that  channel. 

Having  described  how  the  data  is  blocked  and  accessed,  the  overall  flow 
chart  of  IRDSP  can  now  be  outlined  as  in  Figure  4. 

4.  The  Demodulation  Algorithm 

The  radiance  estimation  procedure  used  in  program  TRDSP  is  now  described. 
With  reference  to  the  flow  chart  in  Figure  4,  the  following  description  ap¬ 
plies  to  the  processing  block,  labelled  "ENVELOP." 

Quite  simply,  the  radiance  estimate  is  the  envelope  of  the  complex, 

narrow-band  component  of  the  received  data,  within  the  frequency  band 

f  ±  B  Hz.  Here,  f  is  the  chopper  frequency  and  B  is  the  assumed  radi- 
c  c 

ance-signal  bandwidth.  The  procedure  for  obtaining  the  envelope  entails 
complex  demodulation  and  is  block  diagrammed  in  Figure  5  for  the  ideal 
noise-free  situation. 

The  rest  of  the  section  is  devoted  to  the  justification  for  and  the  limi¬ 
tations  of  complex  demodulation  as  a  viable  procedure  for  estimating  the 
radiance  message  from  chopper-modulated  received  data.  A  noise-free 
analysis  first  establishes  the  validity  of  the  procedure.  Next,  a  suitable 
noise  model  is  included  in  the  analysis;  it  is  shown  that  an  undesirable 
aspect  of  the  in-band  noise  is  to  non-linearly  distort  the  estimate.  This 
affect  will  be  negligible  over  most  of  the  rocket  data  due  to  a  sufficiently 
high  Signal-power  to  in-band  Noise-power  Ratio  (SNR)  and  is  significant 
only  at  very  low  message  levels. 
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Figure  4.  Structured  Flowchart  of  Program  IRDSP. 
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Figure  5.  Block  Diagram  of  Complex  Demodulation  Procedure 
for  Noise-free  Signal. 


Finally,  the  specifications  of  the  practical,  complex  demodulation 
algorithm  are  summarized. 


Noise- Free  Analysis 

The  chopper-modulated  data  channels  are  assumed  to  have  the  following 
form: 


s(t)  =  r(t)  -<p(t)  , 


(7) 


where  r(t)  *  the  radiance  measured  by  the  detector  -  assumed  to  be  low 
frequency  in  nature  and  bandlimited  to  B  Hz 
<£>(t)  =  the  chopping  waveform 


=  2>n  ' 


jnwct 


w  =  2ff  f  =  the  fundamental  radian  frequency  of  the  chopper 
c  c 

[An}  are  the  complex  coefficients  of  the  Fourier  series. 
j6c 

Let  =  A^  •  e  .If  the  received  data  x(t)  is  noise-free, 

x(t)  =  s(t)  =  r(t)  ♦  ^  Ant 
n 

Then,  from  Figure  5,  the  radiance  estimate  is 
r(t)  =  lz(t)  |  =  (  {x(t)  •  e  c  }  *  h(t) 


jnwct 


(8) 


(9) 


(10) 
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Assuming  for  the  moment  that  the  low-pass  filter  is  ideal  and  has  band¬ 
width  jS  >  B,  the  message  bandwidth,  then  one  obtains 


r(t)  «  I  A, 


r(t)  I  = 


Achr(t) 


j9 


=  Achr<4> 


(11) 


i.e. ,  the  estimate  r(t)  is  a  scaled  replica  of  the  radiance  impinging  on  the 
detector.  Note  in  Eq.  (11)  that  the  unknown  phase  of  the  chopper  funda¬ 
mental  is  eliminated  by  the  magnitude  operator,  an  advantage  that  accrues 

jwct 

from  the  use  of  the  complex  sinusoid  e  in  the  frequency  down-transla¬ 
tion  (in  contrast  to  the  use  of  a  real  sinusoid  sin  w  t  or  a  co-sinusoid 

c 

cos  w  t). 
c 


The  Noise  Model  and  Its  Effect  on  the  Estimator 

The  earlier  noise- free  analysis  establishes  the  viability  of  complex  demodu¬ 
lation  as  a  procedure  for  recovering  a  low  frequency,  bandlimited  message 
from  a  noise-free  chopper  modulated  version  of  it.  A  variety  of  noise 
sources,  however,  contaminate  the  actual  modulated  signals.  For  sim¬ 
plicity  these  are  separated  into  two  parts.  One  part  is  the  broadband  ran¬ 
dom  noise  caused  by  background  radiance,  the  electronics  and  the  telem¬ 
etry  transmission  channel.  The  other,  also  assumed  wideband,  is  the 
transient  interference  caused  by  cosmic  ray  showers.  The  received  signal 
can  then  be  expressed  as: 


x(t)  =  s(t)  +M(t)  +w(t)  ,  (12) 

where  s(t)  =  modulated  signal  defined  in  Eq.  (7), 

M(t)  ~  broadband  noise  representing  the  sum  of  background  atmos¬ 
pheric  radiance,  the  detector  noise  and  the  system  noise, 
w(t)  =  cosmic  ray  interference  consisting  of  randomly  occurring  spiky 
transients. 

Note  that  the  two  contaminants  arc  very  different  in  the  time  domain  but 
are  similar  (broadband)  in  the  frequency  domain.  For  this  noisy,  received 
signal,  the  complex  demodulation  procedure  gives 
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(13) 


r(t)  =  j  |x(t)  •  e  ^  c  }  *  h(t) 
j8c 

Ach* 6  •r<t>+Vb.(t)+Wn.b.(t) 

where  an  ideal  filter  has  again  been  assumed  and  n  ,  (t)  and  w  ,  (t)  are 

n.b.  n.b.'  ' 

complex  narrow-band  components  offx(t)  and  w(t),  respectively,  defined 
by  the  spectra 

“n.b.'O  *  u« 

=  0 

and  similarly, 

Wn.b.(f)  =  W<*>  for  fc"^  ^  f  5  fc  +/S 

=  0  elsewhere  .  (15) 

Note  that  both  Un  ^  (f)  and  ^  (f)  are  asymmetric  spectra  since  the  cor¬ 
responding  time  signals  u  ,(t)  and  w  ,  (t)  are  complex. 

n*  u  n.  d* 

On  comparing  Eqs.  (11)  and  (13),  we  see  that  when  the  received  data 
is  noisy,  the  estimate  r(t)  is  not  simply  a  scaled  replica  of  the  radiance 
message,  but  is  corrupted  by  the  in-band  noise  components.  Of  course, 
this  in-band  noise  energy  is  bound  to  show  up  in  the  output,  since  spectrally, 
it  is  indistinguishable  from  the  message  energy.  What  is  undesirable  about 
the  form  of  Eq.  (13),  however,  is  that  the  in-band  noise  distorts  the  mes¬ 
sage  in  a  non-linear  way.  To  see  this  more  clearly  we  proceed  as  follows: 
Let 


for  f  -0  ^  f  <f  +0 
c  c 

elsewhere 


(14) 


j0  (t) 

n(t)=^n.b.<t)+wn.b.<t>=t'<t)‘e  ”  •  (16) 

In  Eq.  (16)  the  noise  representation  is  in  terms  of  the  slowly  varying 
envelope  v(t)  and  phase  0n(t)  of  the  narrow-band,  untranslated,  noise 
process.  Then,  Eq.  (13)  gives 
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r(t)  - 


j0 

Ach  *  r(t)  *  e  C  +  n<t> 

CAcj1r(t)  cos  6c  +  v(t)  cos  en(t)} 

+  j  ^  Achr(t)  8in  6C  +  Wt)  sin  8n(t)}|  , 


(17) 


which  simplifies  to 


r(t)  = 


r^(t)  +  i/*(t)  +  2A^  •  r(t)  •  i/(t>  cos  [©c 


V‘>i 


1/2 


(18) 


Equation  (18)  gives  the  radiance  estimator  when  the  received  data  is  noisy. 
By  setting  t/(t)  =  0  in  Eq.  (18),  the  noise-free  estimate  is  r(t)  =  A^  •  r(t), 
which  is  the  same  as  in  Eq.  (11). 

Of  practical  importance,  however,  is  the  condition  of  relatively  low 
noise,  i.c.,  when 


SNR  =  Power-in  ! 

Power  in  n(t)  **  1  * 


To  see  the  nature  of  the  estimator  at  high  SNR,  rewrite  Eq.  (18)  as 

1 1/2 


(19) 


r(t)  -=  Achr(t) 


“2 'A C8c'6n(*,3+7TJfe: 

oh  Achr  (t) 


(20) 


The  last  term  in  the  expression  is  a  random  process  being  a  function  of  the 
input  noise.  At  time  t,  for  the  random  variable, 


x  =  -v  E)  - 
t  ,  2  2... 

Achr  (t) 

if  the  condition  of  (19)  is  true  it  follows  that 


(21) 


Prob  (xt<  6)  -  1  for  0  <  6  «  1  .  (22) 

Hence,  neglecting  the  last  term  in  (20)  and  retaining  the  first  two  terms  in 
the  power  expansion,  (20)  gives 


14 


tW  -  Achr(t, 


1  +  cos  [8-0  (t)} 

Achr(t)  c  n 

Achr(t)  +  v(t)  cos  (6c  "  ®n(t))  • 


(23) 


Hence,  at  high  SNR,  the  estimate  r(t)  contains  a  scaled  replica  of  the  mes¬ 
sage  plus  output  noise. 

The  conclusion  from  the  above  noise  analysis  is  that  the  complex-de¬ 
modulated  envelope  ?(t)  will  satisfactorily  recover  the  radiance  message 
from  the  IRBS  data. 


Practical  Considerations 

In  the  preceding  discussions,  an  ideal  low-pass  filter  had  been  assumed  to 
simplify  the  discussions.  A  practical  algorithm  cannot,  however,  use  an 
ideal  filter  since  its  non-casual,  infinite  impulse  response  only  allows  the 
computation  of  its  output  with  infinite  delay.  To  make  a  usable  filter,  one 
must,  therefore,  relax  the  frequency  domain  specifications  up  to  the  point 
where  the  performance  is  still  acceptable.  Thus,  the  discontinuous  transi¬ 
tion  from  passband  to  stopband  is  modified  to  provide  an  intermediate  tran¬ 
sition  band  and  rather  than  insisting  on  infinite  stopband  suppression  one 
specifies  a  large  but  finite  attenuation.  Of  course,  the  only  effect  of  using 
such  a  non-ideal  filter  is  that  the  higher  message-translates  in  the  modu¬ 
lated  data  and  the  out-of-band  noise  will  not  be  entirely  eliminated  but 
acceptably  attenuated. 

Once  suitable  specifications  have  been  decided  the  designer  must  choose 
between  the  two  generic  filter  types  -  Finite  Impulse  Response  (FIR)  or  In¬ 
finite  Impulse  Response  (IIR)  -  for  his  needs.  The  attendant  trade-off  is 
that  IER  filters  can  meet  magnitude- response  specifications  with  much  lower 
filter  order  than  can  FIR  filters  but  they  invariably  have  a  non-linear  phase 
response  with  delay  while  FIR  filters  are  easily  constrained  to  have  linear 
phase  and  no  delay.  In  the  present  applicatiou  it  was  decided  that  waveform 
symmetry  and  delay-less  processing  were  uncompromisable  requirements. 
Hence,  an  FIR  linear  phase  filter  is  used. 
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The  specifications  for  the  low-pass  filter  used  in  Figure  5  are  as 
follows : 

1.  FIR.  linear  phase,  real  coefficients 

2.  Passband  (±  .0274  dB)  =  .45  Hz 

3.  Transition  Band  =6.0  Hz 

4.  Minimum  stopband  attenuation  =  50  dR. 

The  filter  was  designed  by  the  windowing  method  using  a  Kaiser  1^  -  sinh 
window.  The  resulting  filter  is  close  to  optimal  (in  the  sense  that  it  is 
"nearly-equiripple").  Its  3  dB  bandwidth  is  2.82  Hz  and  filter  order  is  276. 

When  this  filter  is  used  in  the  demodulation  scheme  on  actual  data,  a 
hitherto  unconsidered  interaction  between  the  sharp  cosmic  ray  transients 
and  the  filter's  long  pi'oeessing  memory  (.49  secs)  comes  to  light.  When  a 
sharply  spiked  cosmic  ray  "hit"  is  processed,  the  output  exhibits  the  familiar 
and  highly  undesirable  ringing  effect.  The  narrower  the  spike  is,  the  more 
this  ringing  is  like  the  impulse  response  of  the  filter.  Obviously,  the  larger 
the  spike,  the  more  pronounced  is  the  ringing.  Also,  if  a  succession  of  hits 
occurs  each  within  .49  secs  of  its  neighbors,  the  radiance  estimate  will  be 
obscured  by  a  long  and  persistent  (and  meaningless)  ringing.  Note  that  the 
above  implied  rate  is  2  hits/sec  which  is  considered  entirely  probable. 

Thus,  we  see  that  our  earlier  conclusion,  that  the  radiance  estimator 
of  Figure  5  was  satisfactory,  was  rather  hasty.  The  earlier  noise  analysis 
was  based  on  steady-stale  frequency  domain  characterization  of  the  con¬ 
taminants,  and  since  both  the  stationary  noise  and  the  cosmic  ray  hits  are 
broadband  phenomena  they  were  lumped  together  into  one  noise  process  in 
Kqs.  (16)  and  (17).  It  is  clear  now  that  the  spikes  have  to  be  removed  and 
that  this  has  to  be  done  in  the  time  domain  where  the  two  contaminants  are 
differentiable. 

The  spike  detection  algorithm  has  to  be  some  variation  of  the  basic  ap¬ 
proach  of  thresholding  the  difference  between  the  signal  and  a  smooth  ver¬ 
sion  of  it.  To  set  a  threshold,  the  local  signal  variance  is  needed.  In  the 
chopped  data,  direct  application  of  such  an  algorithm  will  not  work,  how¬ 
ever,  since  the  variance  values  will  be  dominated  by  the  chopper  energy 
and  not  the  spike  energy.  Thus,  the  chopper  energy  must  be  removed  prior 
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ANALYSIS  OF  THE  DETECTOR  SIGNAL  PROCESSOR 


Initiator:  Mr.  S.  Price 

Problem  No. :  4060 
Project  No. :  7670 

SUMMARY  OF  THE  PROBLEM  AND  PRINCIPAL  RESULTS 
Each  element  in  the  Detector  Array  has  an  associated  Detector  Signal 
Processor  circuit  which  performs  a  logarithmic  compression  of  the  dy¬ 
namic  range  of  the  sensor  measurement.  Our  task  is  to  analyze  the  cir¬ 
cuit  and  to  obtain  a  closed  form  expression  for  the  output-to-input  transfer 
characteristic.  This  characteristic  would  enable  one  to  "de-log"  the 
channel  outputs.  Towards  this  end  we  have  determined: 

(i)  The  "snap-shot"  transfer  function  of  the  logging  circuit  as  a  func¬ 
tion  of  the  output  level, 

(ii)  the  approximate  pass-bandwidth  and  gain  as  a  function  of  the  output 
level , 

(iii)  a  closed-form  expression  for  the  output-to-input  transfer  character¬ 
istic. 

The  transfer  expression  that  we  have  obtained  appears  to  have  an 
adequate  functional  form.  However,  when  it  is  written  in  the  peak-to-peak 
input  versus  peak-to-peak  output  form  to  match  the  calibration  data,  it 
contains  a  dependence  on  the  actual  minimum  and  maximum  output  values 
rather  than  on  just  their  difference.  Because  of  this  a  simpler  expression 
was  heuristically  determined.  The  parameters  in  this  expression  were 
optimized,  via  iterative  regression,  for  each  detector  channel.  The  con¬ 
verged  heuristic  transfer  expression  appears  to  be  quite  satisfactory  in 
explaining  the  calibration  data. 
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PROBLEM  DEFINITION 

The  Detector  Signal  Processor  (DSP)  can  be  represented  by  the  block 
diagram  in  Figure  1,  as  a  cascade  of  four  amplifier  stages. 


Input 

buffer 


r 

"Logging" 

amplifier 

V0(t)  Gain  (5.68) 

- ^  and  Level 

Shift  (2.54) 


d(t) 
_ \ 

Output 

L> 

buffer 

Stage  1  Stage  II  Stage  m  Stage  IV 

Figure  1.  Functional  Block  Diagram  of  the  Detector  Signal  Processor. 


The  end  stages  are  unity  gain,  input  and  output  buffer  stages  respec¬ 
tively,  and  are  used  for  the  usual  obvious  reasons.  Stages  II,  III  are  the 
key  elements  of  the  signal  processor.  In  Stage  II,  the  buffered  sensor  out¬ 
put  is  amplified  with  a  constant  gain  when  |v.  |  <  lmv  and  with  a  loga¬ 
rithmically  decreasing  gain,  with  increasing  |v.  j ,  when  |v^|»lmv.  The 
resulting  bipolar  output  vQ  is  then  amplified  with  a  constant  gain  of  5.68 
by  Stage  III,  and  also  level  shifted  by  2.5  volts  to  produce  a  unipolar  posi¬ 
tive  output  d(t) . 

Our  intent  is  to  obtain  an  output-to-input  transfer  expression  of  the 
form, 


v.=f(d)  ,  (1) 

since  it  will  enable  the  "de-logging"  of  the  channel  outputs,  to  yield  the 
original  sensor  output  v^.  For  this  purpose  it  is  sufficient  to  focus  our 
attention  on  Stage  II,  and  seek  a  transfer  expression, 

Vj  =  g(vQ)  (2) 


because  obtaining  (1)  from  (2)  involves  only  a  trivial  change  of  variables 
defined  by 


Equations  (1)  -  (3)  describe  the  analysis  problem.  In  the  next  section 
we  describe  the  particular  analysis  methodology  that  we  employ. 

ANALYSIS  METHODOLOGY 

Figure  2  below  is  the  circuit  diagram  for  Stage  II. 


Fig.  2.  The  Circuit  Diagram  for  the  Logging  Amplifier 

To  analyse  the  circuit  of  Fig.  2.  we  initially  considered  the  usual  ap¬ 
proach  of  applying  Kirchoff's  current  law  (KCL)  at  the  node  at  the  '-ve' 
input  of  the  operational  amplifier.  The  resulting  differential  equation  is 
non-linear  due  to  the  non-linear  feedback  from  the  logging  transsistors 
Tj  and  T By  linearizing  around  particular  output  values,  one  could 
possibly  obtain  a  solution.  However,  at  this  point  wc  chose  to  use  an  alter¬ 
native  approach  because  we  believe  it  gives  better  insight  into  the  circuit 
action.  This  alternative  approach  for  obtaining  a  transfer  expression  is 
best  explained  in  terms  of  its  three  constituent  steps.  These  are: 
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Step  1:  The  Equivalent  Circuit 

Since  the  logging  transistors  are  acting  as  variable  resistors  which 


change  the  gain  of  the  stage  in  a  non-linear  fashion,  the  first  step  is  to 

replace  T-  and  T„  by  their  incremental  resistances  R  and  R 

1  2  J  v , pnp  v , npn 

which  are  exponentially  decreasing  functions  of  the  output  voltage. 


Conditioned  on  a  particular  value  that  the  output  assumes,  a  linear 
analysis  is  performed  which  yields  a  "snap-shot"  frequency  response. 
Moreover,  the  set  of  values  that  v^  can  assume,  defines  a  family  of  fre¬ 
quency  response  curves.  The  time  varying  function,  vQ(t)  can  then  be 
viewed  as  the  output  of  a  system  whose  frequency  response  is  making  tran¬ 
sitions  from  one  member  of  this  family  of  responses  to  another  with  time. 


For  input  v.(t)  band-limited  to  be  within  the  minimum  bandwidth  that 
the  varying  system  can  ever  assume,  the  system  transfer  function  is  ap- 
proximatable  by  a  frequency-independent  but  out  put -voltage  dependent 


"gain"  (4G(vq)),  i.e., 


Hence, 


v.(t)  G(V 


vi  =  G  (vo> '  vo  =  &<vo>  • 


(4) 


This  is  the  desired  output-to-input  transfer  characteristic,  mentioned 
earlier  as  Eq.  (2). 

In  the  next  section,  the  analysis  outlined  in  Steps  1-3  above,  shall  be 
carried  out. 


CIRCUIT  ANALYSIS 

As  outlined  in  the  previous  section,  the  analysis  will  be  performed  in  three 
steps,  as  follows. 
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Step  1:  The  Equivalent  Circuit 

In  Appendix  Al,  the  incremental  resistance,  Rv,  of  a  pop  transistor 
used  in  the  "transdiode"  logging  configuration  is  derived.  The  result  is, 

“ve/eo 

Rv(VE)  =  V10  <5> 

where 

Rq  =  113.93  MCI  (6) 

EQ  =  72.05  mv  .  (7) 


We  can  define  the  incremental  resistances  of  transistors  T  and  T„  in  Fig. 

1  Ci 

2  in  terms  of  Ry(vE)  of  (5)  as, 


and 


)  = 


<*>  for  vQ  S  0 

W  for  vo>0 


Rv ,  npn^V 


_  ,Rv(iv0l)  for  v0i0 


for  vQ  >  0  . 


(8) 


(9) 


Formulas  (8),  (9)  reflect  the  fact  that  of  the  two  transistors  T^  and  Tg, 
only  one  is  active  depending  upon  the  polarity  of  v^. 

Another  simplification  of  the  circuit  in  Fig.  2  is  that  the  T-network  in 
the  feedback  path  can  be  replaced  by  an  equivalent  resistance, 


Rt  =  2R2 


(10) 


For  R^3  10  KC2  and  =  57.47C2,  formula  (10)  gives, 


Rt  =  1.76  Mft  . 


(11) 


Using  Eqs.  (8),  (9),  and  (11),  the  equivalent  circuit  for  Stage  II  of  Fig. 
2  is  as  shown  below  in  Fig.  3. 
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Fig.  3.  The  Equivalent  Circuit  for  the  "logging"  Stage  of  Fig.  2 


2:  The  Family  of  Frequency  Response  Functions 

The  advantage  of  the  equivalent  representation  in  Fig.  3  of  the  logging 
circuit  is  that  instead  of  analysing  the  non-linear  feedback  explicitly,  we 
can  now  perform  a  linear  analysis  of  the  fixed-elements,  linear  circuit 
that  results  for  each  value  that  the  output  assumes.  Thus,  at  vQ=  g  volts, 
a  "snap-shot"  of  the  varying  frequency  response  is  obtained.  Moreover, 
the  set  of  values  that  the  output  can  assume  defines  a  family  of  frequency 
responses  for  the  circuit.  Useful  simplifications  result  from  this  approach. 
Let 


«b 


=  R_  R  R 

T  v,pnp"  v,npn 


=  the  equivalent  feedback  resistance. 


(12) 


Note:  ||  is  the  parallel  combination  operator.  As  explained  above,  at 
vQ  =  g  a  fixed  value,  Rb  in  (12)  has  a  constant  value.  For  the  circuit  of 
Fig.  3,  we  can  then  define  the  two  impedances. 


Z 


f 


1 

jwC1 


+  R, 


(13) 
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(14) 


zb  =  (^)K 

and  the  circuit  can  be  redrawn  as  in  Fig.  4.  The  notion  of  a  *s nap- shot* 
analysis  can  now  be  elaborated. 


The  impedance  representation  of  Fig.  4  has  been  obtained  at  a  particular 
value  of  vQ=  g.  Still  this  does  not  stop  us  from  obtaining  its  transfer  func¬ 
tion,  since  we  can  always  consider  small  enough  inputs  such  that  the  vari¬ 
ation  in  Zb  can  be  arbitrarily  bounded.  It  is  in  tnis  sense  that  the  transfer 
function  obtained  this  way  is  a  "snap-shot"  of  a  varying  frequency  response. 

For  the  inverting  op-amp  circuit  in  Fig.  4,  the  transfer  function  is. 


H(s) 


^■tRi 

_ sRbCl _ 

(l+sCgl^Ml+sC^j)  ' 


The  magnitude  of  the  frequency  response  can  be  obtained  from  H(s)  as. 


where 


HM(f) 


nl/2 


(f2  +  a2)(f2+b2) 


(15) 
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(16) 


The  important  conclusion  from  (19)  and  (20)  is  that  both  and  HM(fM) 
being  functions  of  b,  are  themselves  functions  of  the  output  voltage. 

From  the  above  "snap-shot"  frequency  response  analysis  of  the  circuit, 
we  find  that  "b"  is  the  critical  parameter  of  the  circuit.  From  (18),  b,  in 
turn,  depends  upon  R^,  the  equivalent  feedback  resistance,  and  it  is  this 
critical  variable  that  we  shall  now  analyse. 

The  equivalent  feedback  resistance, 


Rb'RTllRv,pnpK,np„-RTl|Rv'lvOl>  •  <2l> 

since  only  one  transistor  is  active  depending  upon  the  polarity  of  the  out¬ 
put  Vq.  From  this  expression  one  can  immediately  identify  three  distinctly 
different  regions  of  operation  for  the  circuit  as  follows: 
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Region  I  -  Constant  Gain  Region,  In  this  region,  Rv(vQ)  »  Ibp  ,  i.e., 
the  output  level  is  small  enough  that 


W  *  Ro 


Rn 


Then,  R^s^.RT  and  from  Eqs.  (14)  through  (16)  it  follows  that, 


and 
i.e. , 


But, 


b  =  upper  3  dB  cut-off  =  502.3  Hz 

L.  =/ab  =  44.8  Hz 
M 

HM(fM)  =  174.58 


M'  M' 


II  /C  \  _  S  —  _ *  *-* _ 

a  +  b  C2Rb+C1R1 


a  <  w  <  b  . 


and  since  C2Rb<<  CjR^, 


(22) 


„  (I  ,.£lV5> 

C1R1  Rj  ’ 


(23) 


and  Eq.  (22)  becomes 


HM<W) 


(24) 


This  is  a  significant  result.  It  states  that  for  Vj(t)  bandlimited  to  (a,b)min 
the  "gain"  of  the  bandpass  circuit  is  approximately  independent  of  frequency. 
In  other  words,  the  capacitors,  and  Cg,  have  only  a  marginal  effect  on 
gain.  They  primarily  perform  waveshaping  and  high  frequency  noise  sup¬ 
pression. 

Equation  (24)  gives  us  the  output  voltage  dependent  "gain"  that  we  were 
seeking  in  Eq.  (4).  Thus, 
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(25) 


rQ  R*. 

=  G(v0)  ~  H(w)  ~  -  ~ 


v. 

1 


the  -ve  sign  being  due  to  the  inverting  configuration  of  the  operational  am¬ 
plifier.  Therefore, 


i o  e„ , 


R1  /R1  Rl\ 

'1"‘Rb,V°“  \RT+Rv)V°  ’ 


vo  1  _VEo 

Vi  ~  Gx  "  Gg  V0*  10 


(2G) 


where  Gj  =  F^/Rj  and  Gg  =  Rq/R^ 

Region  II  -  Logarithmic  Gain  Region.  This  region  is  characterized  by 
Rv(v0>  K<  ’  i,e*’  outPut  lev£}l  is  large  enough  that 

-lv0l/E 

Rv*V"V°  °«RT  • 

In  this  case,  Rb~Rv(vQ).  Again  b,  f  and  HM(fM)  can  be  computed  from 
Eqs.  (14)  -  (16)  for  a  fixed  vQ.  The  important  observation,  however,  is  that 


<:  r. 

REGION  II  REGION  I 


and,  therefore,  b  and  f^  can  only  be  larger  than  the  small  signal  values  of 
Region  I.  Similarly,  HM(fM)  can  only  be  smaller  than  in  Region  I. 


Region  III  -  Transition  Region.  This  region,  of  course,  corresponds 
to  those  values  of  the  output  for  which  Rv(vQ)  ~  RT-  The  amplifier  has 
neither  constant  nor  logarithmic  gain  and  hence  the  *crm  "transition  re¬ 
gion." 

The  plot  of  Rv(Vq),  R,p,  and  R^  versus  v^,  below,  shows  the  three 
regions. 
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We  now  have  all  the  necessary  information  to  complete  our  task  of  ob¬ 
taining  a  closed-form  expression  for  the  transfer  characteristic  of  the 
logging  circuit.  This  is  done  next  in  Step  3. 

Step  3.  The  Approximate  Output- Voltage- Dependent  Gain  and  the  Output - 
to-input  Transfer  Expression. 

The  family  of  frequency  responses  found  in  the  previous  section  can 
oe  plotted  for  some  typical  values  of  as  in  Figure  6. 

In  this  plot,  note  that  in  the  frequency  band  (4,502.3)  Hz,  the  responses 
are  flat  to  within  3  dB  of  their  respective  maximum  values. 

This  frequency  range,  of  course,  is  the  pass-band  tor  small  output 
levels.  It  follows  that  if  we  consider  input  signals  which  are  bandlimited  to 
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within  (a,  b)mjn  then  for  such  input  the  frequency  response  can  be  approxi¬ 
mated  by  the  maximum  value,  with  a  known  bounded  error  of  £  3  dB 

and  the  transformation  d(t)  =  5.68  v^t)  +  2.5  can  be  used  to  obtain  Vj  vs.  d 
from  Eq.  (26).  Thus, 

v  «  -  _  (d -_2. 5)  .  -(d~2*5)/5.68E0 

i  5.68G1  5.68G2 


v  =  — ^_5 - 1 - ^  +  — 2J5 — 

i  5.68G1  5.68G1  5.68G2 


2.5/5. 68E  -d/5. 68E 

10  0  •  10 


2.5/5.68E. 


-d/5. 68  E. 


5. 68G. 


•  d  •  10 


v.  -  aQ  +  a^d  +  a210  +  a3  ■  d  •  10 


where  eQ  =  5.68EQ 


a0  5.68G, 


ai  “  "  5.68Gj 

„  c  2.5/5.68E. 

„  _  in  u 

a2  5.68G2  '  10 


2. 5/5. 68E, 


a3  5.68G2  *  10 


Finally,  the  transfer  expression  in  terms  of  peak-to-peak  input  and 
output  values  is, 


r.  =  a,d  +  a„  10 
i,p-p  1  p-p  2  <- 


.  .  "dh/e0  .  “^^O 

+  a3  dn10  "  V  10 


i.e. ,  due  to  the  non-linear  form  of  (27)  the  peak-to-peak  transfer  expres¬ 
sion  requires  a  knowledge  not  only  of  d^  =  d^  -  d^,  but  also  of  the  actual 
minimum  and  maximum  output  levels,  d^  and  d^,  themselves.  This  infor¬ 
mation  has  not  been  recorded  in  the  calibration  data.  Therefore,  to  fit  an 
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expression  to  the  calibration  data,  we  heuristically  fit  a  function  of  the 
form  of  (26) 


v.  =  P~P  +  x, 
i.  P“P  x3  1 


P-P 


10 


d  /x0 
P-P  2 


(29) 


The  constants  x.,  x~,  and  x  were  determined  via  iterative  regression  for 
each  channel.  Results  of  the  fit  for  the  channels  SI  through  S9  are  in¬ 
cluded  as  Figures  7-15. 
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APPENDIX  A1 


The  expression  for  Rv,  the  incremental  variable  resistance  of  the  logging 
transistors  can  be  obtained  as  follows:  For  the  grounded-base  configura¬ 
tion  (known  as  the  "Transdiode"  logging  configuration)  the  modified  Ebers 
and  Moll  equation  for  the  collector  current  is 


I 

c 


o- 


■oV 


E 


I 

c 


aNrES' 


/  qVp/KT  \  /qv  /KT  \  / 

b  /  ~ *cs \  0  -fE'osA 


qv  /M.KT 
,  0  J 


where  v„  =  Emitter-base  voltage 

1 It 

vc  =  collector-base  voltage 

a.T  =  forward  current  transfer  ratio 
N 

OTj  =  reverse  current  transfer  ratio 

M.  >  1  are  the  "uncollected"  current  components  that  flow  through 
the  base  circuit 

q  =  1.60219  X  10~19C;  K  =  1.38062  X  10_23J/°K 
T  =°K  =°c  +  273.15 


Since  v  -  0  due  to  the  virtual  ground  at  the  -ve  terminal  of  the  op-amp 
c 


I 

c 


aNIES* 


qvE/KT 


v 


E 


for 


aNrES 


»  1 


Now  o;N  is  close  to  unity  and  fairly  constant  over  the  entire  range  of  cur¬ 
rents.  Neglecting  the  last  term  (which  for  aN  =  .99  contributes  (1/4)  mv 
offset), 


KT  C  C 

vt.  =  —  X  2.303  X  log  —  =  Eft  log  f —  . 
E  q  ‘ES  U  lES 
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Incremental  Itesistano 


Then, 


'av  E  0  KT 

Rv  =  =  incremental  resistance  at  V£  =  2  3Q3I~  =  ~y  > 

*  r*  r* 


where 


Base- Emitter  Voltage 


R  is  plotted  above  as  a  function  of  the  emitter-base  voltage,  v„.  For 
VE  ^  E0>  t*ie  resistance  falls  by  a  decade  for  every  EQ  increase  in  v^,. 

This  fundamental  property  is  being  used  to  produce  the  logarithmic  com¬ 
pression  in  the  circuit.  Observe  that  for  0  <  v^,  <  Eq.  the  resistance  de¬ 
viates  from  the  above  rule,  [n  fact,  Ry  is  infinite  for  vF  =  0  (the  transistor 
is  cut  off). 

Neglecting  the  "1"  in  the  denominator,  which  is  valid  for  V^/E^  >  1, 


R  =  Rn 
v  0 


10 


■VEo 


This  effect  of  the  approximation  is  that  the  true  R^  (solid  line)  is  being 
approximated  by  the  dotted  line  for  0  <  v^  <  EQ. 

We  now  need  IFg  to  compute  RQ.  In  the  specifications  sheet  provided 
us,  the  base-emitter  junction  has  been  characterized,  as  shown  below,  by 
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the  four  points  shown  on  the  vs.  graph  for  T  =  25°C  and  T=90°C. 


Since  the  circuit  operates  at  T=90°C,  we  are  interested  in  I^g 
This  is  the  intercept  of  the  upper  line.  Let 


T  =  T 


and 


yT  *  log  IF 

T  T  =  t 


Then,  we  seek  to  fit  the  line  =  mxT  +  C.  For  T  =  90°C,  the  scope,  m, 
and  intercept,  c,  are  obtained  using  points  #1  and  #2  in  the  junction  char¬ 
acteristics.  The  resulting  equation  of  the  line  is 


^90°C  =  14*22^9x90°C  ~  9,56122  • 

Hence, 

I  =  io"9-56122  =  o.274f,5  X  10~9 

hb  T  =90°C 

and  substituting  into  formula  for  RQ,  one  obtains, 

RQ  =  1.1393  X  108  =  113.93  M Cl  . 

Also.  Ii  =  l/m  =  70.289  mv,  directly  from  the  slope  'm\ 

Compare  this  value  of  EQ  with  the  value  E^=  72.05.  obtained  earlier 
from  the  physical  equations  describing  the  transistor. 
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A  MODEL  FOR  SATELLITE  CHARGING  AT  HIGH  ALTITUDES 

Initiator:  A.  Rubin 

Problem  No.:  4017 
Project  No.:  7661 

Satellites  in  geosynchronous  orbit  have  been  found  to  charge  up  electrically 
during  magnetic  substorms  leading  to  circuit  upsets  and  malfunction  of 
satellites.  To  understand  the  physical  processes  taking  place  it  is  neces¬ 
sary  to  calculate  the  details  of  the  spacecraft  charging  process. 

A  computer  program  was  developed  to  conduct  numerical  simulations 
of  plasma  interactions  with  a  satellite,  assuming  the  satellite  to  be  repre¬ 
sented  by  an  infinitely  long  cylinder  and  treating  the  plasma  as  discrete 
particles.  The  system  is  assumed  to  be  uniform  along  the  axis  of  the  satel¬ 
lite.  A  number  of  features  of  the  actual  satellite  are  incorporated  into  the 
model  with  as  much  flexibility  as  possible  to  accommodate  a  wide  range  of 
conditions . 

PARTICLE  MOTIONS 

The  particles  are  treated  as  discrete  objects  rather  than  as  fluid  elements 
but  do  not  necessarily  represent  individual  ions  or  electrons.  The  current 
design  associates  with  each  computer  particle  a  cluster  of  identical  elec¬ 
trons  or  ions  so  that  the  several  million  particles  of  the  actual  plasma  are 
represented  by  a  few  thousand  computer  particles.  In  contrast  to  an  earlier 
model  with  spherical  symmetry,  the  parameter  relating  the  number  of  real 
particles  for  each  computer  particle  is  a  constant  for  all  particles. 

The  particle  motions  are  determined  by  the  electric  field  produced  by 
the  satellite  and  by  the  particles  themselves.  The  program  allows  the 
satellite  to  be  charged  to  an  arbitrary  voltage  and  also  allows  the  potential 
in  the  plasma  far  from  the  satellite  to  be  set,  thus  influencing  the  rate  at 
which  particles  emerge  from  the  far  plasma  and  enter  the  region  near  the 
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satellite.  The  electric  field  is  calculated  from  the  electric  potential  ob¬ 
tained  by  solving  Poisson’s  equation  in  a  cylindrical  geometry. 

Thus,  the  equation  to  be  solved  for  the  potential  0  is: 

I  J.  L  *4)  +  A  -  zR 

r  ^r  V  >r)  r2  5Q2  * 

assuming  uniformity  is  the  z -direction.  With  a  coordinate  transformation 
S  =  0nr,  this  equation  becomes: 

^20  +  _  -r2p 

bS2  bS2  f 

which  is  the  form  for  rectangular  coordinates,  using  a  modified  density 

9 

r“p. 

This  transformation  is  effected  by  choosing  the  radial  gridpoints  to 
have  uniform  increments  in  bn  r,  rather  than  as  an  explicit  coordinate 
transformation.  A  standard  Poisson-solving  routine  ^  for  rectangular 
coordinate  systems  is  then  employed  to  find  the  potential,  using  the  modi¬ 
fied  charge  density.  This  approach  also  requires  an  auxiliary  routine  to 
solve  for  the  potential  on  the  satellite  itself,  to  provide  a  necessary  bound¬ 
ary  condition  for  the  Poisson  solver.  The  auxiliary  routine  uses  a  mo¬ 
ment  series  expansion  to  determine  the  satellite  potential  in  terms  of  the 
charge  distribution  and  the  constant  value  assumed  for  the  potential  in  the 
far  plasma  region. 

There  are  also  provisions  for  allowing  a  uniform  magnetic  field  which 
has  the  field  aligned  along  the  cylinder  axis. 

Particles  arc  "lost"  from  the  system  by  striking  the  surface  of  the 
satellite  or  escaping  from  the  region  studied,  and  are  added  to  the  system 
oy  a  variety  of  emission  processes  on  the  satellite,  described  below,  or 
by  entering  into  the  region  studied  from  the  surrounding  plasma. 
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SECONDARY  EMISSION 

Ions  or  electrons  which  strike  the  satellite  are  assumed  to  generate  sec¬ 
ondary  electrons,  taking  into  account  both  the  energy  and  angle  of  inci¬ 
dence  of  the  impacting  particle  to  determine  the  yield  of  secondary  elec¬ 
trons. 

The  equation  used  to  evaluate  the  secondary  emission  for  incident 
electrons  is: 


6  =  C, 


1  -  exp  (-(XQE)n) 


(X0E) 


n-1 


exp  (Cj  [1  -cos  0]) 


where  6  =  number  of  secondary  electrons 
E  =  energy  of  incident  electrons 

0  =  angle  of  incidence  with  respect  to  normal  to  surface. 

The  parameters  CQ,  C  ,  and  XQ  can  be  chosen  to  be  different  over  the 
angular  sectors  of  the  satellite  to  reflect  the  properties  of  different  mate¬ 
rials  with  respect  to  secondary  emission. 

For  incident  ions,  the  equation  used  is; 


aQ/Ef  sec  8 

5  f 


where  a^  and  X^  reflect  properties  of  the  materials  in  each  angular  sector 
of  the  satellite. 

For  either  incident  ions  or  incident  electrons,  the  emitted  secondary 
electrons  are  assumed  to  have  a  Maxwellian  velocity  distribution: 

2 

f(v)  « v  exp 


The  velocity  dispersion  "b"  for  the  emitted  electrons  can  depend  on  the 
species  of  incident  particle,  and  its  value  can  easily  be  changed  at  execu¬ 
tion  time. 

The  direction  of  the  velocity  of  the  emitted  particles  has  a  cosine  dis¬ 
tribution  in  angle  with  respect  to  the  normal. 
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The  corresponding  process  for  ion  emission  (sputtering)  has  not  been 
treated,  as  the  magnitude  of  that  process  is  relatively  small. 


BACKSCATTERING 

Inelastic  collision  of  electrons  with  the  satellite  is  treated  as  a  backseat  - 
tering  process,  again  taking  into  account  the  energy  and  angle  of  incidence 
of  the  impacting  particle  to  determine  the  yield.  The  equation  used  here  is 

/  v  \oose 

*  A(vE>y 

where  )3  -  number  of  backscattered  electrons 
E  =  energy  of  incident  electrons 

9  =  angle  of  incidence  with  respect  to  normal  to  surface. 

The  parameters  AQ  and  BQ  can  likewise  be  set  to  different  values  on  the 

different  angular  sectors  of  the  satellite. 

The  energy  of  the  emitted  particles  is  directly  related  to  the  energy  of 

3  2 

tilt  incident  particles,  according  to  f(v)  «v  /E  ,  for  incident  energy  E, 

2 

provided  m  v  /2  £  E,  and  again  a  cosine  distribution  in  angle  is  assumed. 
As  with  the  sputtering  process,  back  scattering  of  ions  is  neglected. 


PHOTOEMISSION 

The  emission  of  photoelectrons  from  the  satellite  surface  is  treated  in  a 
manner  similar  to  secondary  emission,  taking  account  of  the  incidence 
angle  of  solar  illumination  and  the  photoemissive  coefficients  of  the  differ¬ 
ent  sectors  of  th<-  satellite  to  determine  the  photoemission  yield.  Thus, 
the  relevant  formula  is: 


e  =  kQf(cos  a)  , 

where  (  =  number  of  photoelectrons  per  unit  time  per  unit  surface  area 
f  =  solar  illumination  flux 

a  -  angle  of  solar  illumination  with  respect  to  normal  to  surface. 
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The  solar  flux  is  essentially  zero  or  one.  depending  upon  whether  the  sec¬ 
tor  is  in  shadow  or  not.  while  is  a  photoemission  coefficient  which  can 
be  set  to  different  values  for  different  sectors  of  the  satellite. 

Currently,  the  velocity  spectrum  of  the  emitted  photoelectrons  is 
assumed  to  be  a  Maxwellian  with  the  cosine  angular  dependence  previously 
assumed  for  secondaries. 

PARTICLE  BEAMS 

Corresponding  to  the  particle  beams  on  the  satellite,  there  is  a  routine 
which  emits  particles  from  a  localized  region  of  the  satellite  in  a  beam 
pattern.  The  beam  pattern  is  given  according  to  a  1/0^  cos  0/8^  law,  for 
0^0^  0  and  the  velocity  distribution  of  the  beam  particles  can  be 
chosen  to  be  either  monochromatic  (v  =vQ  for  all  particles)  or  Maxwellian 
^f(v)  <*v2  exp  ^-v2/2v2jj.  There  are  effectively  separate  particle  beams 
for  ions  and  electrons,  and,  in  addition  to  0Q  and  v^,  the  location,  aiming 
direction,  and  current  for  each  beam  can  be  specified. 

An  option  is  available  by  which  the  total  charge  emitted  by  the  beams 
can  be  linked  to  the  total  satellite  charge,  thus  influencing  the  potential  on 
the  satellite. 

ELECTROSTATIC  POTENTIAL  CALCULATIONS  FOR  A  CYLINDRICAL 

MODEL 

Introduction 

Recent  investigations  of  spacecraft  interaction  with  the  surrounding  plasma 
have  led  to  the  development  of  a  two-dimensional  model  for  numerical 
simulations.  This  model  treats  the  spacecraft  as  an  infinitely  long  cylinder 
and  studies  the  motion  of  particles  in  the  cylindrical  annulus  between  the 
spacecraft  and  an  arbitrary  boundary  at  some  distance  into  the  plasma. 

To  represent  the  motion  of  the  particles  properly,  it  is  necessary  to  per¬ 
form  calculations  of  the  electric  field  which  is  established  by  the  distribu¬ 
tion  of  charges  in  the  region  and  by  imposed  potentials  on  the  spacecraft. 
This  section  describes  the  potential  calculation  method  that  has  been  im¬ 
plemented  . 
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Poisson’s  Equation 

For  a  system  with  translational  symmetry  in  the  Z-direction  (along  the 
axis  of  the  spacecraft  cylinder),  Poisson's  equation  is: 


(MKS  units) 

be*  f 


where  the  electrostatic  potential  is  0,  and  p  is  the  charge  density  (per 
unit  area).  The  charges  are  effectively  infinitely  long  rods  in  this  repre¬ 
sentation  and  interact  via  a  logarithmic  potential. 

Converting  the  above  form  of  Poisson's  equation  to  suitable  form  for 
a  discrete  grid  gives: 


where  0^.  -  0(r.,6.),  p.  =  p(r1,6.),  6  =  0.-  6.^,  A}  =  r}+1  -  r.,  and 
^i  "  ri  *-2  ~  ri ' 

For  a  uniform  interval  in  the  radial  grid  points,  this  form  produces  co- 
efiicients  of  0^  ,  ±.  1  j  depend  on  i.  However,  if  A^  is  propor¬ 

tional  to  r.,  these  coefficients  are  constant. 

Thus,  letting  A.  -  ar.,  so  that  D{  =  a(2  +  0f)  r.,  Poisson's  equation 
becomes: 

_  20  +0  1  +.  Vj+i"  2Vi  +  ^J-1  -  Uhl 

a2(2  +a)  L  i+1> j  62  p  ‘ 
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This  formulation  corresponds  to  a  transformation  of  the  radial  coordinate 
according  to  S  =  9m  r,  which  gives 

a20  ,  *20  _  p 

*s2  ae2  f 

2 

where  p  -  r  p  for  the  transformed  density  in  (S,  6)  coordinates. 
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The  charge  density  is  determined  from  the  number  of  charges  within 
a  particular  grid  cell,  but  each  charge  is  treated  as  being  uniformly  dis¬ 
tributed  over  a  cell  in  (S, 0)  coordinates,  rather  than  (r,0)  coordinates. 

o 

This  is  implemented  by  assuming  a  charge  density  distribution  p  =  C../r 

for  each  cell,  normalized  according  to  fr*  r  d6  =  q^. ,  and  then  using  the 

average  of  r2p  over  the  cell  for  the  right  side  of  Poisson's  equation.  Thus, 
2 

r  p.  .is  replaced  by  -q.  ./[6  (l  +a)]  in  the  formula  above. 

*  y  J  *  »  J 

Note  that  different  charges  within  a  cell  do  not  interact,  and  only  the 
net  charge  within  a  cell  contributes  to  the  potential.  Thus,  the  grid  pro¬ 
duces  a  potential  that  is  both  "softer"  (reduced  short-range  interaction) 
and  "smoother"  (spatially  averaged)  than  the  actual  interparticle  inter¬ 
action. 


Boundary  Conditions  for  the  Model 

In  solving  Poisson's  equation,  the  boundary  conditions  must  be  specified 
at  some  radius  in  the  far  plasma  (Rg)  and  over  some  surface  on  or  within 
the  spacecraft.  The  spacecraft  itself  is  modeled  as  a  conducting  cylinder 
completely  surrounded  by  a  dielectric  layer.  This  allows  the  charge  dis¬ 
tribution  on  the  surface  of  the  spacecraft  to  be  that  determined  only  by 
particle  impacts  and  emission  mechanisms,  yet  allows  for  the  specifica¬ 
tion  of  a  unique  fixed  potential  for  the  spacecraft.  Conceptually,  therefore, 
the  two  boundaries  are  concentric  cylinders  with  a  constant  potential  on 
each.  However,  neither  potential  is  necessarily  constant  in  time.  In  par¬ 
ticular,  if  the  spacecraft  potential  is  not  fixed,  it  will  "float"  with  respect 
to  the  outer  boundary  potential,  depending  on  the  charge  distribution  in  the 
intervening  region. 

If  the  boundary  condition  is  left  solely  in  terms  of  the  potential  on  the 
interior  conducting  cylinder  (r  =  Rj),  then  the  routine  for  solving  Poisson’s 
equation  must  also  involve  the  conditions  at  two  interfaces:  the  conductor/ 
dielectric  interface  and  the  dielectric/vacuum  interface.  However,  if  the 
interior  boundary  condition  is  transformed  into  a  specification  of  the  po¬ 
tential  on  the  outer  surface  of  the  spacecraft  (r=  R  ),  then  Poisson's  equa¬ 
tion  need  only  be  solved  over  a  homogeneous  region. 
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This  transformation  can  be  accomplished  by  solving  for  the  surface- 
potential  for  the  case  of  an  isolated  charge,  given  the  potential  at  the 
outer  boundary  and  the  specification  of  the  potential  on  the  interior  con¬ 
ductor  ("floating"  or  fixed  potential).  The  solution  for  an  arbitrary  charge 
distribution  is  then  the  appropriate  superposition  of  solutions  lor  isolated 
charges.  (This  is  essentially  a  Green's  function  approach.) 

For  a  charge  q  at  (r^.0^),  with  0  =  #o  at  r  =  R1  (the  inner  conducting 
cylinder)  and  0  =  ^^at  r  =  R^  (the  far  plasma),  the  potential  at  (R  ,  0)  (on 
the  surface  of  the  spacecraft)  is  given  by: 


*"V9)'alVa2Va3  2ff7* 


cos  n(0  -  0^) 


C  +  — 

11  2ffner 


xMRp/iy 

wher<!  J1  ‘  xfMlip/Kg)  ttV/BpI  ■ 


KK-(Cnl!pn  +  ^H  - 

'  '  p_ 

(R./Rp) 

*2  =  (Rp/Rg)  +MVRp>  =  1  ’  a!  ’ 


-*>  (^/iy  fv  (R./R  ) 

a3  “  X~Hn  (Rp/It^)  +  0r>  (K./R  )  ‘  “a2  ^  ’  and 


”  (Rpn+RB)(Rpn-Rfn)-,<(itp“-R;")(Rpn*Rn 


in  principle,  the  summation  should  no  extended  to  n=°°,  but  on  a  discrete 
grid  it  is  restricted  according  to  the  number  of  angular  sectors  N  (=  2r/<5). 
For  the  case  of  a  floating  potential  on  the  spacecraft,  with  a  total  charge 
Q.  on  the  interior  cylinder, 


Q: 

-  V  2ttx7 


\v  \y|  \V 
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Some  modifications  to  the  above  expressions  are  required  to  accom¬ 
modate  the  convention  of  having  the  charges  distributed  over  each  cell, 
rather  than  as  discrete  charges.  This  is  accomplished  by  a  superposition 
of  the  above  solutions  for  a  charge  density 


H  2 

6  Om  (1  +a)  r 

over  the  region  rc/(l+a)  <  r  £  rQ,  ®  <  ®c  +  T^lis  Pr°duces  a  poten¬ 
tial  #(Rp,0),  with 


where  and  a^  are  as  before, 


bj  =  (R/R) 


*<yV 


+  0r>(R./Rp) 


and 


A  =  R2n 
n  p 


B  =  R2n 
n  P 


(Rr-Rf) 

R2n+R2n) 
P  i  / 

_(rp"+4“) 

(r2*  -  R?n) 

( R2n  -  R2n 
\  P  .1  { 

-  X 

+  x 

(r2"  +  k2"' 

\  P  1 

^R2n+n2n) 

(Rpn  +  Rn) 

(r2"-r2"' 
V  p  > 

^R2n  +  R2n] 

The  condition  for  a  floating  spacecraft  potential  then  becomes 


The  solution  is  now  completely  formulated,  given  a  method  of  solving 
Poisson's  equation  in  rectangular  coordinates.  For  this  particular  appli¬ 
cation,  Hockney 's^  Fourier  Analysis/Cyclic  Reduction  Method  was  em¬ 
ployed.  The  Fourier  transform  routine  for  Hockney's  programs  is  also 
used  *o  obtain  the  superposition  of  solutions  for  setting  the  boundary  po¬ 
tentials  at  r  =  Rp. 
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SPLINE  SMOOTHING  OF  RADIANCE  MEASUREMENTS 


Initiator:  Mr.  S.  Price 

Problem  No.  •  3036 
Project  No.:  7670 

The  task  under  this  problem  number  was  to  develop  one-dimensional  and 
two-dimensional  flexible  knot  least  squares  spline  smoothing  algorithms. 
These  smoothing  techniques  were  applied  to  infrared  radiance  measure¬ 
ments  as  a  function  of  galactic  latitude  and  longitude. 

Least  Squares  Approximations  and  Splines 

This  section  deals  primarily  with  one-  and  two-dimensional  least  squares 
spline  approximations  to  univariate  and  multivariate  data,  but  the  principle 
applies  in  general  to  smoothing  measured  data  by  approximating  the  data 
with  f.‘>  v  ptimal  vector  from  a  designated  vector  space.  In  the  case  of 
cubic  splines,  the  designated  vector  space  consists  of  twice  continuously 
differentiable  functions  which  are  locally  cubic  polynomials.  In  the  case 
of  least  squares  polynomial  approximations  the  designated  vector  space 
consists  of  all  polynomials  of  degree  <;  n. 

I .  One-Dimensional  Cubic  Splines 

A  cubic  spline  approximation  of  a  function  f(t)  defined  on  an  interval  I 
consists  of  local  cubic  polynomial  approximations  to  f(t)  over  subintervals 
of  I  in  such  a  way  that  the  local  approximating  curves  over  contiguous  in¬ 
tervals  join  together  in  a  smooth  fashion. 

Definition:  Let  P  =  [x^.  x^, ...  be  a  partition  of  the  interval  [a,b], 

i.e. ,  a  =  x,  <  x„  <  •  •  •  <  x  =  b. 

'll  n 
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A  function  s(t)  which  is  twice  continuously  differentiable  on  [a,b]  and 

is  a  cubic  polynomial  on  each  sub-interval  (Xj,x.+j),  i<  i<n-l  is  called 

a  cubic  spline.  The  x,  ,x„, . . .  ,x_  are  called  knots  for  the  spline  interval 

i  m  n 

[a,b]. 

Let  S(P,  [a,b])  denote  the  set  of  all  cubic  splines  with  knots  P  *  [x^, 

x„, . . .  ,x  }  over  the  interval  [a,bj.  If  P  and  [a,  b]  are  understood,  then  S 
z  n 

replaces  S(P,  [a,b]). 

There  are  two  basic  uses  of  cubic  splines  -  interpolation  and  smooth¬ 
ing.  Cubic  spline  interpolation  consists  of  finding  a  function  s  in  S  satis¬ 
fying 


s(x.)  =  y .  ,  1  £  i  £  n 


(1) 


for  a  given  set  of  pairs  (x-pYj), . . . ,  (*n,yn).  A  fundamental  result  of 
Spline  Theory  states  that  there  exists  a  piecewise  cubic  function  s  satisfy¬ 
ing  (1)  and,  furthermore,  the  solution  s  is  unique  if  the  additional  con¬ 
straints  s^Xj)  =  c  and  s'fx^)  =  d  are  introduced.  The  problem  of  least 
squares  cubic  spline  approximation  consists  of  finding  a  function  s  in  S 
which  minimizes 


Ntjl-y,]2 


(2) 


for  a  given  set  of  data  values  (t  ,y  ),  1  <,  j  <  M. 

Both  the  interpolation  and  the  smoothing  problem  become  more  mathe¬ 
matically  tractable  by  introducing  vector  space  terminology.  We  observe 
that  the  sum  of  two  cubic  splines  is  again  a  cubic  spline  and  a  scalar  times 
a  cubic  spline  is  a  cubic  spline,  that  is,  S  is  a  vector  space.  The  following 
result  (Prenter,  Splines  and  Variational  Methods,  Wiley-Interscience, 

1975,  p.  80)  is  fundamental  for  the  computation  of  spline  fits. 


Theorem:  The  dimension  of  S  is  n +  2. 


Consequently,  if  c^c^  .  ,c  2 

(2)  consists  of  finding  constants  a  a 

1  £ 

where 


is  a  basis  for  S  then  the  solution  of 

2 

, . . . ,  a  wrhich  minimize  L  e.  , 

’  ’  r>+ 2’  i  ’ 
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(3) 


"V 

'V 

Vi> 

C2<V 

Wi>" 

r\  ' 

(2 

= 
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C2(V 

•••  cn+2«2> 

a2 

• 

_€m_ 

yM_ 
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cn+2(tM* 

• 

_an+2_ 

We  can  rewrite  (3)  as  the  matrix  equation,  y  =  CA  +  e,  and  in  this  form, 
the  desired  least  squares  solution  is  given  by 

A  =  (clc)  1  C*Y  . 
or,  alternatively, 

A  =  C+Y  , 

where  C*  =  (cfC^  C*  is  the  pseudo-inverse  of  the  matrix  C.  Thus,  the 
solution  of  the  least-squares  problem  can  be  reduced  to  selecting  a  com¬ 
putationally  convenient  basis  for  S  and  then  finding  the  pseudo-inverse  of 
the  matrix  C  which  is  generated  by  the  basis  functions.  The  basis  func¬ 
tions  found  to  be  most  convenient  are  those  referred  to  as  B-splines.  In 
the  following,  we  introduce  the  definition  of  splines  and  B-splines  using 
the  approach  of  deBoor.  For  deBoor  the  concept  of  spline  is  more  general 
than  that  introduced  above,  but  the  two  notions  of  spline  coincide  for  twice 
continuously  differentiable  cubic  splines.  Following  deBoor,  a  least 
scjuares  spline  approximation  is  developed  for  splines  of  general  degree 
rather  than  an  algorithm  restricted  to  piecewise  cubic  polynomials. 


Dt  finition:  (deBoor,  p.  108).  Let  t  =  (t.)  be  a  non-decreasing  sequence. 

The  i^-normalized  B-spline  of  order  k  for  the  knot  sequence 
f  is  denoted  B.  ^  ~  and  is  defined  by  the  rule 


B. 


,t<x>  =  Vk 


•  ti>  [ti’ 


•VkJ<- 


X) 


*k-l 


for  all  real  numbers  x. 

Notes:  (1)  for  cubic  splines,  k=4 

(2)  [t.,L  +  J . t.^]  *s  a  <^v*<^ec*  difference  and  may  be  com¬ 

puted  in  the  case  of  distinct  knots  by  the  rule, 
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f(U 


Cti'*i+1  •  •  •  •  W(f>  *  hyws-W—VW1 

. _ «w _ . 
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£(ti+k  V(ti+k“ti+l)’“,ti+k  ti+k-l)- 


(3)  the  definition  of  deBoor  permits  knots  to  coalesce;  each  case 
of  equality  of  adjacent  knots  reduces  the  smoothness  require¬ 
ment  at  that  knot  by  1;  the  cases  treated  herein  deal  only  with 
distinct  knots. 


Remark:  Let  the  knots  for  the  space  Sbex.,x„,x„ . x. .  To  generate 

~  i  l*  »j  n 

a  basis  for  S,  the  above  set  of  knots  is  extended  (in  a  somewhat  arbitrary 
manner)  and  the  normalized  B-splines  generated  for  the  extended  set  of 
knots  yielding  the  desired  basis.  More  concretely,  let  A  be  some  positive 
number.  Define  the  extended  knot  sequence  tj ,  , . . . ,  tn+g  by 


t  =  x  -  3  A  ,  t„  =  x1  -  2A  ,  t  =  x  -  A  ,  t  =  x.  . 


1  £  i  £  n 


t  ,  „ . .  =  x  +  j  •  A  , 
n+3+j  n  J  ’ 


l  ^  j  as  3 


and  let  B.  4  ~(x)  be  the  normalized  B-splines  for  the  extended  knot  set 
F  =  [tj.tg,  ....t^^,}.  Set  B(i,x)  =  Bj  ^  ~(x).  As  a  consequence  of  The¬ 
orem  EX.  1,  deBoor,  p.  113,  we  have  the  following  result. 


Theorem:  B(l,x),  B(2,x), . . .  ,B(n+2,x)  is  i  basis  for  S. 

Note:  The  fact  that  t.  ,t„,  t.,,t  ,.,t  t  .  are  arbitrary  is  discussed  in 

-  1’  2  3  nM  n+5’  n*-n 

deBoor,  p.  114. 
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1. 1  Least  Squares  -  One- Dimensioned  Cubic  Spline 

For  observed  data  (t  ,y(t  (t  ,y(t  ))  and  specified  knots  x  ,x_,...,x  , 

a  i  g  g  l.  n 

the  least  squares  cubic  spline  approximation  models  the  data  as: 


Model:  S(t)  ~  a ^  B(l.t)  +  0t9  B(2,t)  +  •••  +  an+2  B(n+2,t)  +  e,  where 
B(l, B(n+2, . )  are  B-splines.  Incorporation  of  observa¬ 
tion  data  into  the  model  yields. 


yfty 

BU.t^  B^tj)  B(n+2,t1) 

ytt2) 

B(l.t2)  B(2,t2)  ...  B(n+2,t2) 

01 2 

_y(te)_ 

B(l,te)  B(2,te)  •••  B(n+2.te) 

an+2_ 

S  *  Ba  +  (. 


The  least  squares  solution  at  of  S  =  Ba  +  e  is  given  by 


1 . 2  Algorithm  for  Solution  of  Least  Squares  Cubic  Spline  Problem 
(i)  Specify  Knots;  (iv)  Read  observation  data; 

yii)  Extend  Knots;  (v)  Write  S.  B; 

tiii)  Define  basis  functions;  (Vi)  Solve  S  =  Ba  +e 
A  computational  consideration  is  that  for  large  data  sets,  the  compact 
support  of  the  basis  functions  can  be  exploited  in  the  implementation  of 
Step  6. 

Direct  implementation  of  the  definition  given  above  for  the  B-splines 
may  lead  to  numerically  unstable  expressions  involving  the  various  differ¬ 
ence  quotients  as  pointed  out  by  deBoor.  He  suggests  rather  that  the  re¬ 
currence  relations, 


i1  > 

v(xn0 . 


tj***Vi 

otherwise 
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X  -t. 
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i+k-1  i 


f  _ v 

-t,  Bi,k-j(x)  +t^tTBi+i,k-i<x) 


i+k  i+1 


be  exploited  in  the  generation  of  the  B-splines. 

RDP  used  the  above  algorithms  to  generate  a  Fortran  program 
which  does  spline  fits  up  to  order  20.  These  spline  smoothing  programs 
have  been  run  successfully  to  approximate  infrared  radiation  measure¬ 
ments  and  to  generate  calibration  curves  for  infrared  detectors  relating 
measured  voltage  to  measured  radiance. 


II.  Two-Dimensional  Spline  Approximation 

The  object  of  this  section  is  to  describe  a  least  squares  error  spline  ap¬ 
proximation  to  surfaces  Z  =  f(x,y).  A  basic  assumption  is  that  the  surface 
is  defined  over  a  rectangular  region  in  the  xy-plane  and  that  a  rectangular 
grid  of  knots  is  appropriate  for  the  approximation  scheme.  The  algorithm 
employed  assumes  that  the  knots  in  the  x-direction  and  the  y-direction  are 
independent.  If  M^y), . . . ,  Mn(y)  are  the  B-spline  basis  functions  in  the 
y-direction  and  L^x), . . . ,  Ln(x)  is  the  corresponding  basis  in  the  x-direc¬ 
tion  then  the  model  is 

z(x,y)  =TT>ijMi(y)  L.(x)  +  e  .  (l) 

The  basic  scheme  in  the  algorithm  to  be  presented  is  to  reduce  the  two- 
dimensional  spline  separately  to  two  one-dimensional  spline  fits  (after  Call 
and  Judd).  To  describe  the  problem,  let 


xrx2,...,xnl  ,  yry2 . yml 

be  the  knots  in  the  x  and  y  directions,  respectively.  Further,  let 

L1(x).L2(x) . Ln(x)  .  Mj(y),  Mgfy) - >Mm0^ 

be  the  corresponding  B-splines.  If  the  observations  are 
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yry2,-..,ym 

^(Xpyj),  ...,z(xn,yi) 
z(x1,y2) - ,z(xn,y2) 


z(x1,ym),...)z(xn,ym) 


then  the  matrieial  representation  of  (1)  using  the  measured  data  is: 


z(xliyi)  •••  zfx^yj) 

z<xry2)  •••  z(Vy2> 


z(x.,,y  )•••  z(x  y  ) 
'  rJnr  n,Jm' 


Mj(yj)  •••  M^yp 

all  a12  aln 

L1(X1>-*  -  L1(V 

= 

Mj(y2)  •••  Mm(y2) 

a21  a22  a2n 

L2(X1>  •••  L2<V 

• 

M-(y  )•••  M  (y  ) 
l'Jm'  mum' 

aml  am2  amn 

_Ln<xl>  •••  vy 

mbolically. 

Z  =  MAL1  . 

(2) 

The  object  of  the  least  square  fit  is  to  find  the  coefficients  a.,  so  that  the 

i.l 

sum  of  the  squares  of  model  values  minus  the  fit  values  is  a  minimum. 
Algebraically,  the  solution  of  (2)  can  be  effected  in  the  following  manner 


(Continued) 
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{(lA)  M1'/}  =  A1 

r,tf{^MtM]  M1'/} 


(3) 


where  is  the  transpose  of  L.  Using  the  notation  M+  for  the  pseudo-in¬ 


verse  of  M,  i.e. , 

1+  -  (mLm)  Ml 


M 


the  solution  (3)  above  can  be  written  simply  as 

I  U  *  l' 

A 


tv*}']’ 


(4) 


For  the  fact  that  M+Z  yields  the  least  squares  error  solution  to  the  prob¬ 
lem  Z  =  MA  +  f  see,  for  example.  Applied  Regression  Analysis  by  Draper 
and  Smith.  To  generate  the  values 


Z(x.y)  =y^a..M.(y)  L.(x) 

the  following  information  is  needed,  in  addition  to  the  coefficients  a.,,  the 
knots  in  the  x-direction  and  on  the  y-direction  and  a  subroutine  which  gen¬ 
erates  the  basis  functions  L  , . . . ,  L  ,  M^, . . . ,  corresponding  to  the 
given  knots. 

RDP  successfully  implemented  the  above  two-dimensional  smoothing 
algorithm  and  applied  it  to  infrared  measurements  in  the  galactic  plane. 
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STUDY  OF  ELECTRON  DENSITY  PROFILES 


Initiator:  Mr.  J.  Aarons 

Problem  No.:  4083 
Project  No. :  4643 

A  magnetic  tape  from  the  Arecibo  Observatory  containing  high  resolution 
electron  density  measurements  was  received  from  the  initiator.  N(h)  data 
for  two  nights  were  analyzed  to  study  the  variation  of  electron  density  as  a 
function  of  altitude.  A  preliminary  study  of  the  data  showed  the  existence 
of  large  and  small  scale  fluctuations.  To  study  the  large  scale  fluctuations 
the  electron  density  was  determined  at  a  sequence  of  distinct  height  inter¬ 
vals  and  plotted  as  a  function  of  time.  A  spectral  analysis  of  these  time 
series  was  effected  using  the  FFT  technique  to  determine  the  period  of 
these  waves.  For  small  scale  fluctuations,  electron  density  profiles  were 
low  pass  filtered  and  the  spectra  of  the  residuals  analyzed  using  the  Welch 
averaging  technique.  Numerical  integration  of  the  N(h)  profiles  was  per¬ 
formed  to  compute  total  electron  content  for  comparison  with  ground  based 
measurements. 
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SIGNAL  ANALYSIS  OF  ROCKET  DATA 


Initiator:  T.  Conley 

Problem  No.:  4854 
Project  No.:  7660 

This  problem  was  to  develop  a  computerized  system,  including  appropriate 
computer  programs  and  operator  instructions,  to  digitally  filter  (smooth) 
unwanted  noise  from  rocket  measured  radiometer  and  photometer  data  (in 
digital  form),  correct  data  for  rocket  aspect,  and  digitally  differentiate 
smoothed  data  to  obtain  volume  emission  rates.  Current  ICECAP  rocket 
radiometer  measured  data  is  contaminated  by  noise  induced  by  the  telem¬ 
etry  recording  and  digitization  processes  and  which  leads  to  erroneous 
volume  emission  rates  when  differentiated. 

Considerable  data  exists  at  OPR  (approximately  10  rocket  flights) 
which  were  analyzed  by  this  technique. 

Sample  data  was  copied  from  tapes  and  programs  were  developed  to 
perform  initial  processing  on  the  data.  Specifically,  a  noise  filter  and 
differentiation  filter  programs  were  written  and  used  on  the  data.  These 
test  results  were  in  excellent  agreement  with  known  results.  The  net 
noise  suppression  obtained  (noise  filter  minus  differentiation)  was  greater 
than  60  db.  Thus,  development  of  a  complete  system  was  undertaken  as 
outlined . 

PROCESSING  PROCEDURE  FOR  EMISSION  RATE  DETERMINATION 
I.  Establish  data  file. 

Edit  data  (replace  calibration  and  noise  "spike"  points  in  data)  using 

optimal  fitting. 

Convert  to  Brightness  Unit  (Ra). 

Plot  resulting  Ror  data  vs.  time. 

Produce  and  plot  FFT  of  data. 
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Van  Rhijn  Aspect  Correction  (Rq  =  R a/'/ a). 

Plot  resulting  Ra  data  vs .  time. 

Produce  and  plot  FFT  of  data. 

II.  Design  a  noise  filter  and  process  data  through  the  filter  (possibly 
separate  filters  for  precession  and  wideband  noise). 

Plot  filtered  data  (vs.  time). 

III.  Differentiate  output  of  noise  filter  (using  a  differentiation  filter). 

Plot  resulting  output  (vs.  time). 

Convert  to  desired  variables  (emission  rate  vs.  altitude)  and  produce 
plot . 

INTRODUCTION 

The  extraction  of  volume  emission  rate  estimates  from  the  raw  data  is 

difficult  since  a  differentiation  of  the  received  signal  is  required.  Since 

differentiation  is  a  very  noisy  process  (e.g.,  d/dt[Asin  (wt)]  =  wAcos  (wt)) 

the  noise  on  the  received  signal  will  completely  mask  the  desired  output 

after  differentiation.  Thus,  the  noise  level  on  the  signal  must  be  greatly 

reduced  to  allow  us  to  obtain  a  meaningful  estimate  of  the  derivative. 

The  approach  taken  (after  aspect  correction)  was  to  digitally  filter  the 

raw  data  to  "optimally"  pass  the  signal  while  suppressing  the  noise/1' 

Here,  typically  at  least,  80  db  reduction  of  noise  level  was  obtained.  The 

(2) 

resulting  output  signal  was  then  differentiated  using  a  digital  filter.  "  The 
differentiation  enhances  the  noise  by  (at  worst)  20  db.  Thus,  the  overall 
filter  (noise  filter  plus  differentiation)  yields  a  net  processing  gain  in  ex¬ 
cels  of  GO  db.  That  is,  the  S/N  (signal -to-noise  ratio)  out  is  at  least  GO  db 
higher  than  the  input  S/N  ratio. 

I'l!.  TK R  SPECIFICATION 

Net  us  consider  a  more  detailed  discussion  of  the  filters  used.  In  general, 
a  numerical  filter  consists  of  a  set  of  "weights"  which  determine  the 
actual  transfer  function  W(f)  of  the  filter.  (The  design  of  a  numerical  filter 
begins  with  establishing  the  shape  of  the  data  window  in  the  frequency  do¬ 
main  which  will  give  the  desired  effect.)  Having  specified  the  theoretical 
transfer  function,  the  remainder  of  the  problem  consists  of  determining 
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the  weights  in  such  a  way  that  the  actual  transfer  function,  or  fre¬ 
quency  response,  approximates  the  desired  one  as  well  as  possible.  A 
perfect  low  pass  filter,  for  example,  would  leave  unaltered  ail  frequency 
components  from  f  =  0  to  the  desired  cutoff  frequency  f(  and  then  would 
suppress  all  frequencies  greater  than  fT  .  The  response  of  an  actual  nu- 

Lj 

merical  filter  can  only  approximate  this  ideal  behavior,  with  the  accuracy 
of  the  approximation  depending  on  the  values  of  various  design  parameters. 
As  in  the  simple  smoothing  process,  a  numerical  filter  is  applied,  such 

that 


M 

y0(t)  73  Wkx(t  +  kAt)  . 
k=^M 

The  filtering  is  accomplished  by  "sliding"  the  filter  along  the  data,  apply¬ 
ing  it  to  2M  +  1  data  points  to  produce  the  filtered  equivalent  of  the  data 
point  which  has  been  multiplied  by  and  then  moving  each  weight  to  the 
next  point  in  the  series  and  repeating  the  operation.  Repetition  of  the 
process  until  all  the  data  in  a  given  run  have  been  covered  produces  a  se¬ 
ries  of  filtered  data  points  which  defines  the  output  function  y^t).  Within 
the  precision  of  the  filter  these  points  will  trace  out  the  input  function  x(t) 
with  the  unwanted  high  frequency  components  removed  (it  a  low  pass  filter 
is  being  used). 

The  basic  form  of  the  filter  model  is  as  follows.  A  nonrecursive  (or 
transversed)  filter  with  linear  phase  is  described  by  a  transfer  function. 


H(Z)  = 


+  Z 
2 


;owTT 

Evaluating  the  complex  variable  at  Z=eJ 


yields  the  frequency  response 


H(e'2ffF)  -  e 


■j  2  rrF  M 


M 


k=0 


cos  27TFk  =  e 


-j2TTFM 


»0<F) 


where  F  =fT  is  the  normalized  frequency  and  T  is  the  sampling  interval. 

With  the  filter  coefficient,  W^,  constrained  to  be  real,  the  frequency 
response  is  composed  of  a  linear-phase  term,  e  and  the  purely 
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real  frequency  response  H  (F)  of  a  zero -phase  filter.  Since  the  linear- 
phase  term  only  introduces  a  finite  delay  of  M  samples,  the  design  prob¬ 
lem  becomes  one  of  fitting  HQ(F),  a  mirror  image  polynomial,  to  the 
desired  frequency  responses.  That  is,  determine  the  coefficients  W 

k 

which  yield  the  desired  frequency  responses  for 

M 

Hq(F)  =y>kcos  2ffFk  . 
k=0 

For  the  case  at  hand  the  weights  were  chosen  to  minimize  the  mean  squared 
error*1)  between  the  actual  filter  frequency  response  and  an  ideal  low  pass 
filter  with  a  cutoff  frequency  chosen  in  accord  with  the  input  signal.  Here, 
a  correction  factor  for  the  "Gibbs"  phenomena  was  included  to  provide  a 
smooth  transition  from  the  pass  band  to  the  stop  band  of  the  filter. 

The  filters  used  here  are  specified  in  the  figure  below.  It  was  found 
that  narrower  filters  began  to  remove  some  of  the  signal  information  while 
wider  filters  allowed  precession  and  more  noise  through  with  no  added 
signal  information.  Thus,  the  bandwidth  was  chosen  to  lie  between  these 
two  conditions. 
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Let  us  consider  the  differentiation  technique  employed.  The  ideal 
frequency  response  characteristics  of  a  digital  differentiator  are  shown 
below.  The  curves  show  the  magnitude  and  phase  of  the  frequency  re¬ 
sponse.  The  magnitude  response  increases  linearly  up  to  a  normalized 
frequency  of  1.0  (the  Nyquist  frequency)  and  then  decreases  linearly  back 
to  0.0  at  the  sampling  frequency.  The  magnitude  response  is  periodic  in 
frequency,  as  shown,  because  of  the  discreteness  property.  The  phase  is 
tr/2  radians  for  frequencies  up  to  the  Nyquist  frequency  and  -it/ 2  radians 
from  the  Nyquist  frequency  to  the  sampling  frequency,  and  is  also  periodic. 

Magnitude 


0  12 


Frequency  Response  Curves  for  an  Ideal  Differentiator 

The  technique  employed  to  approximate  an  ideal  differentiation  filter  is 
given  in  Reference  2. 

In  this  technique,  linearly  spaced  samples  of  the  frequency  response 
of  the  desired  filter  were  specified  and  the  continuous  frequency  response 
was  determined  using  the  discrete  Fourier  transform.  The  interpolation 
formula  obtained  was 


Phase 

It 


2 

0 

1 

2 

I 

ff 

2 


i.e.,  [H^]  are  values  of  the  continuous  frequency  response  at  equally  spaced 
points  around  the  unit  circle;  T  is  the  sampling  period;  and  N  is  the  duration 
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of  the  impulse  response  in  samples.  By  making  the  substitution 


H, 


jffk/N 


(3) 


each  of  the  terms  inside  the  summation  of  (1)  becomes  imaginary,  and  thus 
the  entire  sum  is  imaginary.  For  N  even,  the  complex  factor  outside  the 
summation  in  (1)  represents  a  pure  delay  of  an  integer  (e-jw^N^)  plus 
one-half  number  of  samples.  Thus,  (1)  suggests  that  a  differ¬ 

entiator  with  exactly  half  a  sample  delay  can  be  designed  nonrecursively  by 
setting 


jk/(N/2)  ,  k  =  0, 1, . . . ,  N/2 

v (n  -k)/(N/2)  ,  k  =  (N/2)  +  1, . . . ,  N-l 


and  applying  the  substitution  of  (3)  into  (1). 

Further  details  of  the  design  technique  and  resulting  response  curves 
can  be  found  in  Reference  2. 

For  the  work  described  here  a  differentiation  of  lengths  128  was  em¬ 
ployed. 


SAMPLING  RATE 

When  experimental  data  are  derived  by  discretely  sampling  some  phen¬ 
omenon  at  equally  spaced  intervals  of  time,  the  problem  of  aliasing  may 
occur  in  which  the  sampling  rate  is  low  enough  to  confuse  two  or  more 
frequencies  in  the  data. 

The  net  result  is  that  they  appear  to  be  the  same  frequency.  To  avoid 
this  problem  and  hence  to  define  a  unique  input  function  as  described  by  a 
set  of  data  points,  one  must  be  able  to  assume  that  the  phenomenon  studied 
is  spectrally  limited  to  the  range  |f|<f  ,  where  f  =  fr/2,  f_  being  the 
the  cut-off  or  Nyquist  frequency.  If  such  an  assumption  is  valid,  then  the 
function  has  been  sampled  frequently  enough  so  that  all  significant  fre¬ 
quency  components  are  determinable.  This  is  a  result  of  the  sampling 
theorem  of  information  theory.  The  sampling  theorem  states  that  if  a  func¬ 
tion  G(t)  contains  no  frequencies  higher  than  W  cycles  per  second,  then  it  is 
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completely  determined  by  giving  its  ordinates  at  a  series  of  points  spaced 
1/2W  seconds  apart,  the  series  extending  throughout  the  entire  time  do¬ 
main. 

Since  the  bandwidth  of  the  last  analog  filter  before  digitizing  is  approxi¬ 
mately  100  Hz,  only  about  200  to  400  samples  per  second  are  required  (al¬ 
lowing  for  the  analog  filter  roll  off). 

Furthermore,  since  the  desired  signal  bandwidth  was  confined  to  less 
that  0.1  Hz  (as  can  be  seen  from  the  FFT  plot),  considerable  aliasing  can 
be  tolerated.  That  is,  sampling  rates  of  30  to  50  samples  per  second  will 
only  cause  aliasing  above  0. 1  Hz  thus  not  affecting  the  signal  but  simply  in¬ 
creasing  the  noise  in  the  higher  frequency  range.  With  60  db  processing 
again  this  noise  increase  can  be  tolerated  allowing  a  lower  sampling  rate. 

In  fact,  a  sampling  rate  of  33  samples  per  second  was  used  for  the  proces¬ 
sing  described  here. 
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WAVELENGTH  PEAK  STUDY 


Initiator:  T.  Conley 

Problem  No. :  4855 
Project  No. :  7670 

This  problem  was  to  develop  a  computerized  system  to  smooth  noise  con¬ 
taminated  spectral  features  in  rocket  measured  CVF  data  (ICECAP  data) 
in  order  to  accurately  specify  wavelength  of  peak  emission  features  and 
determine  line  widths. 

A  noise  filter  was  developed  for  this  data.  Also,  techniques  for  "op¬ 
timal"  estimation  of  peak  locations  were  being  studied.  In  particular,  the 
possible  use  of  full  band  differentiations  (with  variance  estimates)  was 
studied. 

Due  to  high  noise  levels  and  short  data  spans,  it  was  found  that  the  re¬ 
quired  accuracy  could  not  be  obtained  by  these  techniques.  It  was  decided 
that  further  analyses  was  not  warranted  at  this  time.  Hr  '  ever,  the  noise 
filter  developed  was  made  available  to  the  initiator. 
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STATISTICAL  ANALYSIS  OF  GAS  MEASUREMENT  DATA 


Initiator:  C.C.  Gallagher 

Problem  No. :  4861 
Project  No.:  6687 

The  problem  was  to  develop  statistical  techniques  to  identify  gas  measure¬ 
ments  from  column  and  spectrometer  tests.  Using  calibration  and  test 
runs,  procedures  for  specifying  an  identification  procedure  was  to  be  de¬ 
veloped  including  probability  measures  of  the  uncertainties. 

A  program  was  developed  to  calculate  the  "statistics"  of  the  calibra¬ 
tion  and  test  runs  and  provide  student  t  tests  for  most  probable  compounds. 
Sample  data  was  used  to  evaluate  the  approach.  The  program  is  described 
below. 

STRATOSPHERE  GAS  MEASUREMENT  ANALYSIS  PROGRAM 


1.  Column  Measurement 

A.  Establish  calibration  reference  base  for  time. 

B.  Process  experimental  data  (time  of  arrival  and  concentration). 

(i)  Take  as  many  runs  as  possible. 

(ii)  Establish  statistical  distribution  of  time  and  concentration 
(e.g.,  averaging  independent  samples  of  time  of  arrival  yields 
a  Gaussian  estimate  with  a  variance  dependent  upon  the  data 
stability). 

C.  Compare  experimental  with  calibration  data  to  establish  "most 
likely"  compounds.  Here,  one  can  order  possible  compounds  ac¬ 
cording  to  their  probabilities.  [Note  that  closest  in  time  is  not 
necessarily  most  probable.] 

2.  Process  spectrometer  data  in  similar  fashion. 

3.  Check  if  the  results  of  1  and  2  agree  as  to  most  probable. 
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d)  If  not:  perhaps  need  more  data. 

A  second  program  was  developed  to  plot  and  fit  gas  concentration 
measurements.  Specifically,  gases  captured  during  a  balloon  flight  were 
chemically  analyzed.  The  data  from  the  analyzer  was  handpunched  ac¬ 
cording  to  a  certain  format. 

This  program  plotted  gas  concentration  versus  amplitude  at  different 
pressures,  and  then  made  separate  least  square  fits  at  each  pressure. 
The  plots,  therefore,  aided  the  user  in  extrapolating  characteristics  of 
the  gas  and  of  the  equipment  at  undetectably  low  pressures. 
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SIGNAL  STATISTICS  OF  SCINTILLATION 


Initiator:  H.E.  Whitney 

Problem  No. :  4893 
Project  No.:  4643 

For  this  problem  RDP,  Inc.  functioned  as  mathematical  and  statistical 
consultants.  Here,  we  answered  a  number  of  questions  concerning  the 
"best"  procedures  for  processing  and  analyzing  scintillation  data. 
Following  is  a  summary  of  the  important  results. 

COMMENTS  AN1)  QUESTIONS  CONCERNING  THE  MAJOR  ASPECTS  OF 
THE  DATA  ANALYSIS 

We  shall  present  some  comments  and  questions  concerning  the  data  anal¬ 
ysis  techniques  which  should  be  employed  for  this  problem. 

Before  considering  specifics  we  make  the  following  general  remarks: 

(1)  When  the  sample  size  of  a  data  stream  is  large  it  is  not  always  neces¬ 
sary  (or  advisable)  to  use  "optimal"  processing  procedures.  In  fact, 
using  "reasonable"  approaches  often  yield  results  very  close  to  op¬ 
timal  while  providing  simplicity.  However,  in  general,  one  should 
know  how  one  is  performing  relative  to  "optimal"  so  as  to  justify  the 
procedures  employed. 

For  the  problem  at  hand  where  the  system  bandwidth  is  2Hz  and 
the  data  samples  are  15  minutes  long  one  has  approximately  3600  in  ¬ 
dependent  samples  per  data  record.  In  general,  this  provides  a  large 
number  of  degrees  of  freedom,  allowing  the  use  of  suboptimal  tech¬ 
niques. 

(2)  The  desired  outputs  are  critical  in  determining  the  techniques  to  be 
employed.  In  particular,  the  required  resolution  bandwidth  and  con¬ 
fidence  bars  are  key  factors. 


I .  Estimate  of  Auto  Correlation  Function 

The  appropriate  estimator  for  the  auto  correlation  function  is  given  by: 
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This  provides  an  unbiased  estimator  with  a  known  variance, 
cally.  the  variance  is  approximated  by: 
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(T  - r-  |u  I ) (R*  (u)  +  R(u  -  r)  R(u  +  r))  . 
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The  exact  statistics  (p.d.f.)  of  Rf  are,  in  general,  not  known.  However, 
the  variance  provides  a  measure  of  the  quality  of  the  estimate. 

Should  M  be  large,  one  must  consider  obtaining  Rf  by  employing  the 
FI  T  approach  (Faster  and  better  estimator  for  large  M).  Here,  one  must 
consider  adding  zeros  to  the  data  stream  (N  zeros  needed  for  aliasing)  and 
employing  weighting  to  provide  a  stable  estimate. 


if.  Rower  Spectra  Estimate 

The  power  spectra  estimate  may  be  obtained  as  (a)  the  Fourier  transform 

of  the  auto  corr  -.1  at  ion  function  or  (b)  directly  by  use  of  the  FFT  on  the 

original  data  stream.  The  major  considerations  here  are  as  follows: 

(a)  Fourier  transform  of  R 
- — - -  r 

(i)  Desired  rc  hition  Bandwidths  (13^)  and  Degrees  of  Freedom 
(DOF) 

With  no  further  weighting,  we  have 
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where  AT  =  time  between  samples. 

Further  weighting  leads  to  the  following  relations: 
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These  different  windows  have  varying  characteristics  which 
do  not  appear  to  be  critical  for  our  problem.  What  is  im¬ 
portant  is  to  obtain  the  desired  B  and  DOF. 

e 

(ii)  Statistics  of  Spectral  Estimator 

o 

The  resulting  spectral  estimator  is  \  distributed  with  the 
DOF  given  above.  Thus,  confidence  bars  can  be  established 
accurately. 

(b)  Spectra  by  Direct  FFT 

(i)  Desired  resolution  Bandwidths  (Be)  and  Degrees  of  Freedom 
(DOF) 

Using  the  direct  FFT  without  any  time  and/or  frequency 
averaging  leads  to  a  poor  spectral  estimator.  In  fact,  the 
standard  error  is  equal  to  the  estimate  itself.  It  is  not  con¬ 
sidered  critical  how  the  averaging  is  done.  That  is,  the 
primary  factors  are  the  resulting  bandwidth  resolution  and 
the  DOF  obtained.  Thus,  frequency/time  averaging  or  use 
of  some  weighting  function  is  somewhat  arbitrary  and  we  can 
consider  the  options.  It  is,  however,  important  to  correctly 
determine  the  DOF  which  result  since  this  dictates  the  ap¬ 
propriate  Chi-squared  statistics. 
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General  Comments 


(i)  Bias  vs.  variance  (fidelity  vs.  stability) 

Smoothing  of  the  estimator  in  order  to  reduce  the  variance 

causes  a  bias  error.  Essentially,  smoothing  increases  the 

resolution  bandwidth  B  which  reduces  the  variance  while  in- 

e 

creasing  the  bias  error.  The  DOF  obtained  essentially 
establishes  the  variance,  while  the  bias  error  is  given  by: 

Ki 

Bias  =  — r*(f)  , 

2 

M 

where  r*(f)  =  second  derivative  of  the  spectra 

K  K2  =  constants  depending  on  the  smoothing  procedure. 
Some  general  properties  of  the  bias  error  are  as  follows: 

•  When  (f)  is  negative  (near  a  peak  of  the  spectra),  the  bias 
is  negative  and  thus  peaks  will  be  underestimated.  Con¬ 
versely,  when  T'"(f)  is  positive  (near  a  trough),  the  bias 

is  positive  and  thus  troughs  will  be  overestimated. 

•  The  narrower  the  peaks  or  troughs,  the  larger  the  value 
of  Tvr (f)  and  thus  the  larger  the  bias. 

•  The  bias  is  reduced  as  M  increases  (the  base  width  of  the 
window  is  decreased).  Therefore,  the  wider  windows 
have  larger  bias  error  (and  smaller  variance). 

(ii)  Correlation  of  samples 

Vriien  smoothing  is  employed  the  samples  become  correlated 
(resolution  bandwidth  increases!.  One  should  realize  that  the 
adjacent  samples  are  not  independent.  However,  the  samples 
should  be  plotted  since  it  is  possible  to  miss  a  peak  whose 
frequency  lies  halfway  between  uncorrelated  samples  if  only 
the  independent  samples  arc  plotted. 

II i .  Estimates  of  Cumulative  Probability  Amplitude  Distribution  (C.D.F.) 

If  it  is  appropriate  to  plot  the  empirical  cdf,  the  procedure  for  the  present 
problem  should  be  t.o  group  the  data  (i.e.,  set  up  a  frequency  table)  and 
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plot  the  cumulative  distribution  for  the  sample.  There  is  too  much  data  to 
plot  the  empirical  cdf  directly.  Problems  that  arise  here  are  as  follows: 

(a)  determining  how  best  to  group  the  data  so  that  we  will  not  distort 
the  data  (e.g.,  should  groups  be  equi -probable  or  fixed  width  with 
a  standard  set  of  class  limits  ?) 

(b)  determining  if  there  are  additional  ways  to  display  the  data  so  that 
it  is  informative  (e.g.,  histograms). 

Also  in  presenting  the  empirical  cdf  is  it  useful  to  put  con¬ 
fidence  limits  around  it  so  that  a  graph,  such  as 


can  be  presented. 

(c)  Goodness  of  fit  tests. 

It  appears  that  there  are  two  main  families  of  distributions  to  con¬ 
sider— the  Nakagami  m-distributions  and  the  log  normal  distri¬ 
butions.  At  present,  the  test  used  for  goodness  of  fit  is  the  chi- 
squared  test.  The  points  that  must  be  considered  here  are: 

(i)  determining  the  number  of  categories  to  use  in  the  chi -squared 
test  that  will  give  the  test  maximum  power. 

(ii)  determining  how  best  to  estimate  the  parameters  of  the  dis¬ 
tribution  under  consideration.  In  order  to  use  the  chi-squared 
test  all  the  unknown  parameters  of  the  distribution  being  con¬ 
sidered  must  be  estimated  from  the  data.  There  are  optimal 
ways  to  do  this  (e.g. ,  maximum  likelihood  techniques).  Due 
to  the  large  sample  size  sub-optimaJ  techniques  (e.g., 
methods  of  moments)  may  be  sufficient. 

(iii)  determining  the  appropriate  degrees  ci  freedom  for  the  ehi- 
squared  test. 


(iv)  One  should  also  consider  alternative  testing  procedures.  For 
example,  the  Kolmogorov  Smirnov  test  might  he  more  pow¬ 
erful  than  the  chi-squared  test.  Also,  the  and  b9  tests 
applied  to  test  for  the  log  normal  are  more  powerful  than 
either  the  chi-squared  or  Kolmogorov  Smirnov  tests. 

(v)  There  may  be  situations  where  more  than  one  distribution 
may  fit  the  data.  For  example,  the  Nakagami  m-distribution 
with  m  =  1  may  fit  the  data  while  the  Nakagami  with  m  =  2  may 
also  supply  a  good  fit.  We  may  be  able  to  supply  techniques 
for  deciding  which  of  these  two  is  the  better  fit  of  the  data. 

In  addition,  we  may  be  able  to  supply  good  graphical  tech¬ 
niques  to  aid  in  deciding  which  distribution  best  fits  the  data. 

IV.  Other  Considerations 

(a)  High  frequency  power  spectra  slope  estimate 

The  estimate  of  the  slope  of  the  high  frequency  roll  off  of  the 
spectra  should  be  studied. 

Most  probably  a  least  square  fit  to  the  spectra  will  provide 

"best"  results.  However,  we  should  consider  other  alternatives. 

2 

For  example,  the  statistics  of  the  spectra  are  known  (X  dis¬ 
tributed)  and  the  slope  is  really  an  estimate  of  the  derivative  of 
the  spectra.  Thus,  we  may  be  able  to  establish  the  statistics  of 
the  slope  and  from  this  the  maximum  likelihood  estimate  (which 
may  be  approximately  or  exactly  the  least  square  fit). 

(b)  Aliasing 

Since  the  system  bandwidth  is  2Hz  and  the  sampling  rate  is  6  sam¬ 
ples/sec,  there  should  be  no  real  aliasing  problems.  When  this 
is  a  concern,  it  can  be  checked  by: 

(i)  Increasing  the  sampling  rate,  or 

(ii)  Checking  the  high  frequency  level  of  the  spectra.  This  value 
should  reflect  the  thermal  noise  level  of  the  system.  (Look¬ 
ing  at  some  samples  this  appeared  to  be  satisfied.) 
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(c)  DC  level 

The  DC  level  should  be  removed  from  the  data  before  further 
processing  is  done.  This  can  cause  numerical  differences  in  es¬ 
timates  (for  example,  the  bias)  if  not  removed. 

(d)  Maximum  entropy  spectral  estimates 

When  interested  in  estimating  a  power  spectra  with  sharp  peaks, 
the  maximum  entropy  spectral  estimate  has  pi'oven  to  be  superior 
to  the  Fourier  estimate.  The  spreading  and  bias  problems  in¬ 
here  in  the  Fourier  approach  limit  its  use  for  peaked  spectra. 

CHI-SQUARED  TEST  FOR  GOODNESS  OF  FIT 

The  chi-squared  test  can  be  used  to  test  if  a  set  of  data  is  a  random  sam¬ 
ple  from  a  specified  distribution.  We  suggest  using  this  test  with  equi- 
probable  cells.  The  number  of  cells  should  be  from  about  10  to  20.  The 
final  number  wall  have  to  be  decided  after  discussion. 

The  steps  in  the  test  are  as  follows: 

(1)  Decide  upon  k,  the  number  of  cells. 

(2)  Select  the  hypothesized  distribution,  F(x). 

(3)  Find  the  values  of  x},  such  that  if  the  true  distribution  is  F(x), 

F(xj+1)  -  F(x.)  =  i  .  xQ  =  -» 

Xk  =  *  * 

These  x.  determine  the  k  cells. 

i 

(4)  Compute  the  number  of  observations  in  the  cells,  call  these  0. 
for  i  =  1, . . .  ,k. 

(5)  Compute  the  chi-squared  statistic, 

— 9 

(6)  If  F(x)  is  the  time  distribution,  X“  has  a  chi-squared  distribution  with 
k-1  degrees  of  freedom. 
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Often  times  there  are  parameters  of  F(x)  which  are  not  known  and 
must  be  estimated  from  the  data.  Say  there  are  m  of  these  unknown  pa¬ 
rameters.  These  parameters  should  be  estimated  using  the  maximum 
likelihood  estimation  procedure.  All  of  the  steps  (1)  to  (5)  above  are  car¬ 
ried  out  as  shown  except  F(x)  of  (2)  has  the  estimates  in  it  and  not  the 

—2 

actual  parameters.  The  distribution  of  X  of  (5)  is  now  chi-squared  with 
k-m-1  degrees  of  freedom  when  F(x)  is  the  true  distribution. 

The  IMSL  (International  Mathematics  and  Statistics  Library)  subrou¬ 
tines  can  be  used  to  perform  some  of  the  above. 

Subroutine  GFIF  can  be  used  to  compute  the  chi-squared  statistic. 

The  hypothesized  F(x)  must  be  given  as  input.  All  parameters  (or  esti¬ 
mates  of  them)  must  be  supplied.  The  number  of  estimates  must  be  ap¬ 
plied  as  input. 


LSK  OF  KOLMOGOROV  -SMIRNOV  TEST  TO  OBTAIN  CONFIDENCE  IN- 
TKKVAL  FOR  EMPIRICAL  CUMULATIVE  DISTRIBUTION  FUNCTION 
It  is  possible  to  predict  by  means  of  confidence  intervals  how  close  the 
cumulative  distribution  of  a  sample  can  be  expected  to  be  to  the  cumula¬ 
tive  distribution  of  the  population.  The  confidence  interval  is  based  on 
the  Kolmogorov-Smirnov  test. 

The  steps  in  the  procedure  are  as  follows: 

(1)  Arrange  the  data  in  order  of  magnitude 


(2)  Plot  empirical  cumulative  distribution  function 
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For  our  problem  we  must  investigate  if  the  data  must  be  grouped. 
Possibly  with  the  rounding  in  the  data  grouping  will  not  be  needed. 

(3)  Above  and  below  the  empirical  distribution  draw  two  parallel  polygons 
at  a  distance  100(1. 36/n)  =  136/n  (for  a  .95  confidence  interval).  The 
width  of  a  band  giving  confidence  .95  for  the  statement,  "The  cumu¬ 
lative  distribution  of  the  population  is  in  the  band"  is  200(1. 3f>/n)  = 
272/n. 

(4)  The  confidence  band  will  look  like 


Any  population  cumulative  distribution  function  which  fits  completely 
in  this  shaded  region  is  consistent  with  the  data  at  the  .95  confidence 
level  (or  .05  significance  level). 

The  IMSL  (International  Mathematics  and  Statistics  Library)  subrou¬ 
tines  can  be  used  to  perform  the  above. 

Subroutine  VSORTA  can  be  used  to  sort  observations. 

Subroutine  USPC  can  be  used  to  print  and/or  plot  empirical  cumulative 
distribution. 

Subroutine  USPC  can  be  used  to  print  and/or  plot  confidence  interval 
(either  a  .95  or  .99  confidence  intervals  can  be  obtained). 

Subroutine  USPC  can  be  used  to  print  and/or  plot  a  theoretical  (popu¬ 
lation)  empirical  cumulative.  This  theoretical  distribution  can  be  plotted 
on  the  same  graph  with  the  empirical  distribution  function  and  the  confi¬ 
dence  band.  It  is  possible  to  produce  printouts  and/or  plots  for  as  many 
theoretical  distributions  as  the  user  desires  (e.g.,  log-normal.  Nakagami 
m-distributions). 

There  are  a  few  problems  that  may  occur  with  the  above  procedure. 

(1)  The  sample  size,  n,  may  be  very  large.  This  may  create  computa¬ 
tional  problems. 


(2)  The  Kolmogorov -Smirnov  test  and  the  resulting  confidence  interval 
assumes  a  continuous  distribution  function  and  no  grouping  or  rounding 
in  the  sample  data.  The  technique  is  not  valid  if  these  conditions  are 
not  met. 

(3)  The  theoretical  distributions  that  can  be  appropriately  compared  to  the 
empirical  distribution  (and  the  confidence  band)  are  supposed  to  be 
completely  specified.  There  are  not  supposed  to  be  any  estimated 
parameters.  The  technique  is  not  valid  if  the  parameters  are  esti¬ 
mated  . 

(4i  It  may  be  possible  that  no  mathematical  model  (theoretical  distribution) 
really  matches  the  data.  If  the  sample  size  is  large  we  may  observe 
the  phenomenon  that  no  standard  distribution  (e.g.,  log-normal  or 
Nakagami-m)  is  acceptable  by  the  above  statistical  criteria  (i.e, ,  falls 
within  the  confidence  band). 

None  of  the  above  problems  should  be  taken  to  imply  that  this  approach 
should  be  discarded.  What  is  implied  is  that  a  routine  writing  of  computer 
programs  will  not  produce  a  finished  product.  Further  discussion  is 
needed  to  decide  how  to  handle  these  problems.  Also,  it  must  be  realized 
that  questions  will  arise  as  the  programs  art  being  put  together.  These 
questions  will  have  to  be  answered  as  they  arise. 

MAXIMUM  LIKELIHOOD  ESTIMATES  OF  THE  PARAMETERS  OF  T1IE 

NAKAGAMI-m  DISTRIBUTION 

One  way  of  writing  the  Nakagami-m  distribution  is 

f(s)  -  —  M—  Sm-1c*ms/fi 

Omr<m) 

where  S  =  signal  power.  0  -  average  power. 

<  ®  (m  =  1  for  Rayleigh  fading,  m  =  l/2  for  fading  more 
severe  than  Rayleigh),  and 
T(m)  -  gamma  function  of  m. 


H4 


Here, 


f2  =  f  sf(s)  ds  ,  m  = 
J0 


/( sf(s)  ds)2  1 

f(s  -Cl)2  f(s)  ds  2 


The  likelihood  for  a  random  sample  of  size  n  is 

,.mn  -(m/0)ES. 

L=— ^ - fS  S„-  •  *S  ,m_1  e 

nmnrn(m)  12  " 


a  log  l 
3m 


=  n  log  m  +  nm  —  -  n  log  fl  -  n  log  S.  — — - 

*»  r(m)  6  i  n 


3  log  L  _  -nm  _m 
aO  Q 


-nm  m  V~^ 

TT  tfl/i 


The  solution  to  these  equations  after  setting  3  log  L/3m  =  5  log  L/3fi  =  0 


%  es  rr/,m>  i  sioga 

ft  =  ___  =  s  and  [^-logmJ  =  ~ir~-logS  . 
The  last  equation  must  be  solved  recursively. 

N.B.  r'(m)  . 


INTERPRETATION  OF  SCINTILLATION  INTENSITY  SPECTRUM 
Representing  the  scintillation  random  process  as  R(t),  one  may  use  either 
inphase  and  quadrature  or  magnitude  and  phase  to  describe  the  process. 
That  is.  y(t)  | _ 


x(t)  ~  j  R(t)  I  cos  6( t) 
y(t)  =  t H(t)  |  sin  0(t) 
R(t)  =  I  R(t)  |  ej6(t) 


6(t)  i 

- 1 - x(t) 
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Since  only  the  scintillation  intensity  is  measured  (i.e.,  0(t)  unknown),  our 
signal  is  ( H(t)  |2.  Classically,  the  power  spectrum  is  interpreted  as 

P(w>  =  J[0(w)]  =  !E[.?[R(t)]}  |2  =  |Rp(w)|2  , 

(watts  If  it  exists  for 

or  db  watts)  each  realization 

of  the  random  process 

where  0(w)  is  the  statistical  correlative  function  of  R(t)  and  E[  }  denotes 
the  mathematical  expectation.  However,  here  we  define  the  spectra  of  the 
scintillation  intensity  as 


Pg(w)  =J[|R(t)i2] 


for  a  given  realization  of  the  process. 

Using  the  relation  of  multiplication  and  convolution  of  the  Fourier  trans¬ 
forms,  we  have 


Pg(w)  =J[iR(t)|2]  -  J[R(t)  R*(t)] 
=  Rp(w)  f  R*(w)  =  $g(w)  , 


!  ition  process. 


Thus,  P  (\v)  is  the  spectra  of  the  scintillation  intensity  (in  units  of  (db 

watts))  and  is  not  the  same  as  the  power  spectra  P(w)  (or  the  square  of  the 
2 

power  spectra  P<w)  )  of  the  scintillation  process  itself.  But  P  (w)  is  re- 

s 

luted  to  P(w).  In  fact,  it  is  the  correlation  function  of  the  spectra  Rp(w). 

P  (w)  can  be  interpreted  as  a  spectra  by  considering  the  scintillation  in¬ 
tensity  as  the  time  signal  of  interest. 

.,2  2 

Finally,  |Ps(w)  1  would  have  units  (db  watts)  and  would  be  the  power 

spectra  of  the  scintillation  intensity. 
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SPECTRAL  SHAPE  ESTIMATE 

Attempts  tu  obtain  the  MLE  (Maximum  Likelihood  Estimate)  for  the  spectral 

slope  lead  to  an  intractable  mathematical  formula  due  to  the  fact  that  the 

2 

spectral  estimates  are  Chi-squared  (\  )  distributed  rather  than  Gaussian. 
Thus,  we  recommend  use  of  the  least  squares  estimate  (which  is  the  MLE 
under  the  Gaussian  assumption).  Here,  the  possible  models  are  shown 
below. 


Case  (a)  requires  a  simple  LS  fit.  However,  cases  (b)  and  (c)  require  a 
choice  of  w.,  the  break  frequency.  Here,  we  suggest  choosing  w.  to  mini¬ 
mize  the  overall  mean  squared  error  over  w.  Reference  (1)  provides  a 
procedure  for  determining  w. .  Basically,  the  procedure  involves  entering 
our  initial  guess  for  w.  and  evaluating  equation  (1)  below  in  a  neighborhood 
about  this  point.  Here,  the  final  choice  for  w.  is  that  value  which  maxi¬ 
mizes  L(i). 


L(i)  =  N  log/2?  -  i  log  CTj  -  (N -i)  log  a2- ^  ,  (1) 

where  N  =  total  number  of  frequency  points  where  spectra  is  estimated  and 
ct  ,  a  L  the  standard  deviations  of  the  LS  estimates  of  the  spectra 

JL  Li 

for  the  sections  below  w.  and  above  w  respetively.  for  the  i  Ixnng  evalu¬ 
ated. 
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FADE  DURATION  DISTRIBUTION 

To  determine  the  distribution  of  fade  duration,  define  Fn(t)  by 


#(t<t  > 

Fn(t)  = - — .  ,  n  ~  #  fades  of  a  specific  magnitude. 

Wo  can  then  assume  the  fades  are  independent  separate  samples  and 
ihe  fade  distribution  is  an  empirical  distribution  function. 

Since  the  true  (theoretical)  distribution  function  is  not  known,  we  can 
only  establish  confidence  intervals  by  using  a  non-parametric  procedure. 
We  suggest  the  Kolmogorov-Smirnov  procedure.  Here,  the  confidence  in¬ 
tervals  are 


Fn<V 


to 


Fn<V+d 


where  (if  the  sample  size  n  is  greater  than 50),  we  have 


For 

confidence 

d 

0.  HO 

1.07//n 

0.50 

1.22//S 

0.95 

l.;J6//n 

0.99 

1 . 63//n 
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INVERSION  OF  CONVOLVED  III  DATA 


Initiator:  Thomas  D.  Conley 

Problem  No. :  4910 
Project  No. :  7070 

The  major  effort  done  under  this  problem  number  was  the  simulation  and 
analysis  of  two  optical  (IR)  spectra  detection  systems.  One  involved  an 
electronic  Butterworth  filter  while  the  second  replaced  the  Butterworth 
configuration  with  a  temporal  filter.  The  equations  may  be  summarized 
as  follows: 

1.  Using  Butterworth 

System  Transfer  Function  is  given  by 


UNCLASSIFIED 
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2.  tTs 


Filter 


Here,  the  Butterworth  filter  is  replaced  by 


Performance  of  these  systems  with  a  target  signal  applied  is  evaluated 
by  replacing  B(f)  with  the  target  spectra. 

The  relative  performances  of  these  two  systems  may  be  evaluated  by 
comparing  the  output  S/N  ratios  of  the  two  systems. 

The  parameters  are  defined  as  follows: 


X  4  footprint  in  meters 
Vj  4  drift  velocity  m/sec 
N  ^  order  of  Butterworth  filters 
Y’t  4  target  velocity 
M  4  1.3 
K  4  constant 


=  V  /2X 

=  fj/10 
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ALIGNMENT  EXPERIMENT  DATA  ANALYSIS 


Initiator:  Capt.  J.A.  Shearer 

Problem  No. :  4921 
Project  No.:  7628 

This  problem  was  best  handled  using  multiple  regression  analysis  for  fitting 
theoretical  curves  to  experimental  data.  Consequently,  the  quality  of  various 
physical  models  were  assessed. 

In  this  particular  application,  the  microscopic  tilt  of  a  missile  in  a  silo 
was  used  as  the  dependent  variable,  as  a  function  of  fifteen  other  independent 
variables  such  as  time,  temperature,  other  tilts,  and  their  derivatives. 

The  program  calculates  the  correlation  among  the  independent  vari¬ 
ables,  and  correlation  with  the  dependent  variable  and  the  residuals  result¬ 
ing  from  the  fit  are  analyzed  to  give  their  auto-correlations  and  power 
spectra.  This  allows  analysis  of  any  structure  that  may  remain  if  the  model 
is  incomplete. 

Results  of  the  analysis  showed  that  «  90%  of  the  variation  in  micro¬ 
scopic  tilt  was  accounted  for  by  an  approximate  function  of  temperature 
and  certain  tilt  and  tilt  rate  measurements. 
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PICTURE  ELEMENT  WORD  LENGTH  REDUCTION  STUDY 


Initiator:  Rupert  S.  Hawkins 

Problem  No.:  4929 
Project  No.:  8628 

INTRODUCTION 

The  problem  concerned  the  development  of  a  program  for  obtaining  efficient 
computer  oriented  techniques  for  compressing  the  word  length  of  satellite 
imagery  while  retaining  Image  integrity.  What  we  had  to  start  with  was  an 
array  of  image  brightness  values,  G(x,y).  The  array  is  rectangular,  and 
each  of  the  brightness  values  in  the  array  can  be  associated  with  a  small 
rectangular  area  of  the  image  called  picture  elements. 

The  purpose  of  this  effort  was  to  obtain  efficient  automatic  computer 
techniques  for  transferring  image  information,  spread  over  the  conven¬ 
tional  six  or  eight  bits  of  an  image,  to  one  bit  while  retaining  the  same 
element  array.  The  objectives  were  the  preservation  of  image  brightness 
information  in  the  reduced  version  and  fast  execution  of  the  techniques  on 
conventional  computers. 

The  above  program  takes  into  consideration  that  the  original  image, 
G(x,y),  can  be  either  six  or  eight  bits  per  picture  element.  Results  were 
obtained  for  the  following  array  sizes:  1)  80x80.  2)  500x500,  and  31 
6400x6400  elements. 


Mathematical  Formulation  and  the  Optimal  Solution 

.he  objective,  in  the  ideal  case,  for  this  problem  involves  solving: 


mm  z 


ZEE 

x  y  n=l 


G(x,y)n  -  ^(x,y)n  P 


where  G(x,y)n  =  average  of  values  at  n  grid  points  about  (x.y)  using  6  or  8 
hits  information  (64  or  256  levels) 
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and  ^(x^yj^  =  is  to  be  chosen  so  as  to  minimize  z  as  a  similar  average 
using  1  bit  of  information  (2  levels). 

The  problem  as  stated  is  an  integer  programming  problem  with  a 
quadratic  objective  function/**  The  sizes  of  standard  arrays  (grids)  to  be 
considered  are  80  x  80 ,  500x500,  6400x6400.  Note  that  this  indicates,  in 
the  smallest  case,  that  6400  variables  are  to  be  assigned  integer  values. 
Even  in  this  "small’case  we  are  well  into  the  realm  of  "large"  mathe¬ 
matical  programming  problems.  In  the  absence  of  the  integrality  condition 
the  problem  is  a  nonlinear  program.  On  the  surface,  it  may  appear  the 
added  requirement  of  integer  valued  variables  should  not  present  a  serious 
problem.  Indeed,  the  solution  space  is  more  restricted  and  the  search 
need  not  be  made  over  an  infinite  number  of  points.  Unfortunately,  such  a 
conclusion  could  not  be  any  more  erroneous.  Although  the  solution  space 
for  the  integer  (discrete)  problem  is  structurally  better  defined,  it  has 
proved  to  be  computationally  formidable.  In  spite  of  over  two  decades  of 
intense  research,  and  a  tremendous  increase  in  computer  speed  and  power, 
the  developed  integer  algorithms  have  not  yielded  satisfactory  computational 
results. 

Some  algorithms  have  been  developed  for  nonlinear  objective  integer 
programming  problems,  which  are  only  of  theoretical  interest  at  the  cur¬ 
rent  time.  Even  for  small  problems,  18-20  variables,  they  are  not  satis¬ 
factory.  For  large  problems,  over  100  variables,  the  state-of-the-art  is 
severely  limited  even  in  theoretical  results,  computational  advances  are 
virtually  non-existent. 

For  these  reasons,  this  optimization  problem  cannot  be  solved  with 
the  present  state-of-the-art  in  integer  programming.  However,  we  report 
this  approach  formally  here  since  it  defines  the  mathematical  structure  and 
offers  some  hope  for  the  future. 

Additional  discussion  of  integer  programming  is  provided  in  Appendix 
A  for  the  sake  of  completeness. 


General  Considerations  From  the  Literature 

A  review  of  the  literature  in  picture  processing  and  pattern  recognition* 
has  led  to  the  following  basic  conclusions: 
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•  There  appears  to  be  no  other  work  closely  related  to  the  idea  of  con¬ 
serving  local  brightness  as  a  criteria  for  bit  reduction. 

•  Virtually  all  algoritnms  for  picture  processing  are  somewhat  heuristic. 
That  is,  they  do  something  reasonable  and  consider  the  effects  of  vary¬ 
ing  the  algorithms  rather  than  formally  deriving  an  "optimal"  approach. 
This  is  due  largely  to  the  size  of  the  arrays  generally  of  interest  and 

to  the  inherent  nonlinearity  of  the  desired  optimization. 

•  Generally,  considering  the  problem  of  thresholding  or  truncation  and 

(2) 

"enhancement,"  the  basic  forms  of  the  algorithms'  '  involve  the  use  of 
the  Laplacian  with  smoothing  and  subsequent  thresholding.  Here,  de¬ 
pending  on  the  particular  objective,  the  details  of  a  specific  algorithm 
vary. 

It  is  important  to  note  that  the  present  algorithms  agree  with  this 
general  approach  and  are  thus  consistent  with  the  literature. 

Based  on  these  findings  and  the  difficulties  associated  with  integer  pro¬ 
gramming,  it  is  felt  that  the  basic  form  of  the  algorithm  is  proper.  Thus, 
the  effort  should  be  (and  has  been)  directed  toward  the  computer  logic  and 
parameter  value  optimizations  of  the  algorithm,  rather  than  attempting  to 
develop  a  "new*  approach. 

Computer  Optimization  Efforts 

We  present  here  a  discussion  of  the  changes  considered.  We  include  both: 

1.  the  historical  or  developmental  approach  of  what  has  been  done,  and 

2.  the  conclusions  and  recommendations  resulting  from  the  effort. 

Historical 

Briefly  stated,  an  initial  examination  was  made  of  the  program,  and  a  pre¬ 
liminary  List  of  possible  places  for  improvement  was  made.  A  few  changes 
were  made  immediately  to  reduce  the  core-memory  size  of  the  program, 
even  before  any  bench-marks  were  run.  Then,  a  considerable  amount  of 


¥ 

Including:  NASA  Reports  and  Defense  Documentation  Center  Reports. 
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time  was  spent  in  getting  SPY  operating  properly.  (SPY  is  a  timing  program; 
it  differs  from  SECOND  in  both  accuracy  and  focus  of  attention;  the  program 
being  timed  itself  is  essentially  unchanged. )  Next,  our  own  timers  and 
counters  were  inserted  in  the  program.  Lastly,  subroutine  MATH  was  tested 
under  a  variety  of  frame  sizes,  shapes,  and  relative  positions,  with  a  variety 
of  changes  in  Fortran,  machine  language  and  compiler  options. 

Conclusions  of  Computer  Code  Optimization 
For  all  picture  sizes: 

1.  Subroutine  MATH  was  improved  by  about  40%  by  changing  the  structure 
of  all  "IF"  statements  logic  and  using  computer  optimizations  level 
OPT=2.  Further  optimization  to  the  "unsafe'1  level  shows  no  essential 
improvement. 

Examination  of  the  machine  codes  shows  that  rewriting  the  present 
logic  in  machine  language  will  give  no  further  improvement. 

2.  In  practical  applications  on  the  CDC  6600,  subroutine  MATH  is  much 
less  important  than  previously  believed  (when  considering  the  overall 
program  running  time). 

For  small  pictures  (80  x80  elements)  almost  all  the  time  is  spent 
in  PACK  and  ROW;  even  if  we  were  able  to  reduce  the  time  spent  in 
MATH  to  zero,  it  would  have  little  effect  on  the  overall  running  time. 

For  large  pictures  (672x475  mcidas)  still  twice  as  much  time  is  spent 
in  PACK  as  in  the  improved  MATH. 

Under  the  present  setup,  frames  which  are  wide  and  short  run  faster 
than  tall,  narrow  ones,  because  the  number  of  calls  to  PACK  is  directly 
proportional  to  the  height.  Once  in  PACK  the  entire  horizontal  line  is 
processed,  not  matter  how  little  is  actually  needed.  We  recommend  trying 
new  ways  of  logic  in  PACK,  such  as  random  access,  or  keeping  pointers 
to  the  starting  edge  of  the  frame  desired.  A  destred  frame  size  of 
6400X6400  looks  as  if  it  might  cause  overload  problems  when  considering 
the  practical  capabilities  of  the  CDC  (>600  system  (by  "system"  we  mean 
the  conglomerate  of  hardware,  software,  and  computer  operations  as 
practiced  at  AFGI.).  A  few  things  to  be  examined  further  are: 


(a)  fine  tuning  of  buffer  sizes 

(b)  find  out  why  PACK  and  ROW  are  better  when  not  optimized 

(c)  program  should  be  changed  to  give  a  message  if  the  frame  size  is 
bad,  instead  of  running  to  time  limit  with  no  printout. 

(d)  the  only  things  on  the  original  checklist  that  remain  to  be  con¬ 
sidered  are  buffering  and  double  buffering.  These  should  be 
attacked  after  the  scheme  for  PACK  is  decided  on. 

The  Two  Algorithms  Resulting  from  the  Study 
The  two  algorithms  will  be  referred  to  as: 

1.  Original  Optimized. 

2.  Modified. 

The  original  optimized  algorithm  is  that  due  to  R.  S.  Hawkins  which  we 
have  optimized  for  fastest  computer  running  time  while  yielding  the  iden¬ 
tical  final  matrix  (0,1).  A  flow  diagram  is  shown  in  Figure  1. 

The  Modified  algorithm  changes  the  logic  yielding  a  different  final  re¬ 
sult.  This  algorithm  is,  in  general,  faster;  however,  in  some  cases  the 
"quality"  of  the  final  (0, 1)  picture  is  inferior  to  that  of  the  original  optimized 
algorithm.  A  flow  diagram  for  the  modified  algorithm  is  given  in  Figure  2. 
The  major  difference  is  the  skipping  of  elements  which  have  been  previously 
pushed  to  the  upper  or  lower  bound. 

General  Conclusions 

We  found  that  the  results  of  both  algorithms  were  very  dependent  upon  the 
values  of  the  incremental  step  size  (c  in  Figures  1  and  2)  and  the  upper 
and  lower  levels  (a  and  b  in  Figures  1  and  2). 

The  Modified  algorithm  had.  in  general,  Improved  running  time  and 
'•educed  energy  difference.  However,  the  "picture  quality"  seems  inferior 
in  some  cases  (for  example,  more  granular). 
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*  GO  TO  DIAGONAL  OR  IF  IN  DIAGONAL  MOVE  NEXT  POINT 

Fig.  1.  Original  Optimized  Algorithm 
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APPENDIX  A 

Some  Properties  of  Integer  Programming;  Problems^ 1  * 

Standard  Integer  Programming  Problem  with  linear  objective, 

n 

max  or  min  >  c  x 

l-j  i  J 

j=l 

m 

subject  to  ^~^a..x.  =  b.  ,  i  =  l,2,...,m 

J--1 

x.  ^0  ,  j  =  l,2,....n 

x.  an  integer  (for  some  or  all  i) 

J 

1.  Cannot,  as  a  first  step,  ignore  the  integer  value  constraints,  solve  the 
resulting  math  progi-am  (linear  or  nonlinear)  and  then  round  off  the 
answers.  Reasons  are  (1)  rounded  values  may  not  satisfy  the  con¬ 
straints  (i.e..  not  feasible),  or  (2)  rounded  value  if  not  necessarily 
even  close  to  the  optimal  integer  solution. 

Cannot  solve  continuous  problem  and  then  round  up/down  the  answers. 

2.  Methods  of  solution  generally  fall  into  one  of  two  approaches: 

i)  enumeration  techniques  including  branch  and  bound  techniques  or 
ii)  cutting  plane  techniques. 

That  is. 

Branch  and  Bound  Method  -  Idea  is  to  partition  the  feasible  region  into 
more  manageable  subregions  -on  each  one  of  these  branches  a  bound 
on  the  objective  is  derived  and  when  it  becomes  clear  that  the  Branch 
will  not  lead  to  the  optimal  solution  attention  is  shifted  to  another  branch. 

There  are  a  number  of  explicit  and  implicit  enumeration  techniques 
which  have  been  generated.  Such  techniques  are  generally  in  the  class 
"branch  and  bound  algorithms." 

Note  that  if  there  are  n  variables  which  are  allowed  to  assume 
only  two  values  each,  then  a  complete  tree  has  2n  nodes  (in  the  picture 
processing  "small  case",  we  have  2^^  in  the  80x80  case,  2*^  in 
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10x10  case).  As  a  consequence  if  n  is  large,  the  amount  of  computer 
time  required  is  very  large.  In  large,  binary  cases  an  implicit 
scheme  is  used.  Further,  in  linear  objective  case,  at  each  branch 
point  a  linear  program  (continuous  variables)  must  be  solved. 

Cutting  Plane  Methods  —  are  almost  always  outperformed  by  a  branch 
and  bound  method,  ft  works  by  modifying  linear  programming  solu¬ 
tions  until  the  integer  solution  is  obtained.  It  does  not  partition  the 
feasible  region;  rather  it  uses  a  single  linear  program  which  is  modi¬ 
fied  by  adding  new  constraints.  It  is  a  very  inefficient  algorithm. 

Each  new  constraint  (cut)  removes  an  L.P.  non-integer  solution  region 
and  when  a  solution  is  reached  which  contains  only  integral  values  it 
must  be  the  optimal  point. 
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COMPUTATION  OK  OFF-AXIS  RADIATION 

Initiator:  Mr.  B.  Schurin 

Problem  No.  :  49-19 
Project  No. :  7670 

A  computer  program,  OARRJ,  was  written  to  calculate  the  off-axis  radia¬ 
tion  from  the  earth's  atmosphere  measured  by  a  sensor  as  a  function  of 

(i)  sensor  altitude 

(ii)  tangent  height,  and 

( i i i >  wavelength  interval. 

The  program  assumes  that  LOWTRAN  4  has  been  used  in  emission 
mode  to  compute  atmospheric  radiance  with 

(a)  a  fixed  sensor  altitude,  and 

ib)  twenty- five  lines  of  sight  covering  zenith  angles  of  90’’  to  180°. 

Also  assumed  is  that  the  spectral  intervals  and  their  corresponding 
filter  functions  have  been  specified  and  that  these  filter  functions  have  been 
convolved  with  the  output  radiance  from  LOWTRAN  4.  Currently,  this  com¬ 
putation  is  performed  in  a  program  called  LOWFILT. 

input  to  OARFJ  consists  of: 

(a)  LOW  KILT  output  spectral  radiance  incident  on  the  sensor  for  25 
lines  of  sight,  and 

|b)  sensor  off-axis  rejection  weights. 

The  product  of  the  sensor  off-axis  rejection  function  and  the  spectral 
radiance  of  the  earth's  atmosphere  is  numerically  integrated  over  a  solid 
angle  in  OARFJ.  Since  the  off-axis  rejection  function  is  computed  relative 
to  the  sensor  line  of  sight,  the  coordinate  system  used  in  LOWTRAN  to 
compute  the  radiance  is  rotated  to  yield  radiance  values  compatible  with 
the  new  coordinate  system  before  the  integration  is  performed  in  OAREJ. 

A  subroutine  is  included  in  OAREJ  to  accommodate  different  vehicle  alti¬ 
tudes.  This  option  permits  the  user  to  run  LOWTRAN  at  a  fixed  vehicle 
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altitude  and  then  obtain  off-axis  results  for  a  variety  of  different  vehicle 
altitudes.  Alternatively,  LOWTRAN  may  be  run  for  each  vehicle  altitude 
and  then  LOWFILT  and  OAREJ  run  using  each  of  the  LOWTRAN  outputs. 

Description  of  the  Computations  in  OAREJ 

In  Figure  1,  the  unprimed  axes  define 
the  sensor  coordinate  system  and  the 
primed  axis  define  the  coordinate  sys¬ 
tem  in  which  the  radiance  from  the 
earth's  atmosphere  is  calculated.  The 
origin  of  the  two-coordinate  systems 
is  the  sensor  altitude  (vehicle  altitude); 
the  //-axis  is  the  local  zenith  and  the 
z-axis  is  the  sensor  line  of  sight;  fi  is 
the  angle  from  the  local  zenith  to  the 
Figure  1  line  of  sight. 

The  transformation  from  the  sensor  to  the  radiance  coordinates  can 
be  written  as 


/ x'\  /  cos  )?  0 
I  y'  1*1  0  1 

\z'J  \-sin)3  0 


(1) 


In  terms  of  spherical  coordinates  (p,  0,<p),  where  p  is  the  distance,  0  the 
azimuth  angle,  and  <0  the  zenith  angle,  the  transformation  can  be  rewritten 
as 


psinip'  cos  0'\  /  cos  )9  0  sin  0 

psin  <p'  sin  0'  |  0  1  0 

poos  <p'  y\-sinj3  0  cos)? 


psin  <p  cos  0' 
psin  sin  0 
pc  os  <o 


(2) 


Only  the  last  equation  in  (2) 
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cos  <p'  =  cos  0  cos  0  -  sin  0  sin  0  cos  0 


(3) 


is  required  in  the  off-axis  radiation  calculations  since  the  radiance  in  the 
atmosphere  is  assumed  to  be  independent  of  the  azimuth  angle. 

In  the  sensor  coordinates,  the  off-axis  radiation  measured  by  the  sen¬ 
sor  as  a  function  of  0,  the  angle  from  the  zenith  to  the  line  of  sight,  is 
given  by 


OAR(  0)  =  j  dfi  W <0,  0)  I<0'(0 , 0  ,0))  , 


where  dfi  =  sin  <n  ckn  d0  =  differential  solid  angle; 

W(0,0>  =  point  source  off-axis  rejection  function  of  the  sensor 

2 

1(0')  =  radiance  of  the  earth's  atmosphere  (watts/cm  -ster) 
1(0')  is  a  function  of  <n,8  and  0  as  defined  in  Eq.  (3),  i.e. , 


i«>,e,0>  . 


Equation  pi)  can  be  rewritten  as 

f  w/2  f  ff/2 

OAR( 0)  =  I  d0sin0  I  W(0 , 0)  [(0,0,0)  d0  .  (6) 

J(D ,  ■'-ff/2 

d 

in  the  Kq.  ;(i),  the  integration  over  0  in  the  present  calculations  has 
been  restricted  to  the  lower  hemisphere  about  the  line  of  sight.  Numerical 

integration  of  Eq.  (6)  is  performed  in 
0  o  7‘  the  computer  program  using  the  trap- 

ezoidal  method. 

As  shown  in  Figure  2,  the  zenith 

^  angle  0  and  the  tangent  height  TH  are 

T»\  related  by 

\  RE  TH+R 

E\  sin  (IT  -  0)  =  sin  0  =  —  ,  (*) 

\  E 

Figure  2 
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where  H  =  the  altitude  of  the  sensor,  and 
R£  =  the  radius  of  the  earth. 

In  the  computation  of  OAR(/0)  in  Eq.  (f»)  the  first  problem  is  to  deter¬ 
mine  the  values  of  I(<p,0,j8).  The  known  values  are  i(<p')  computed  using 
LOWTRAN  4  and  LOWFILT.  The  relation  between  o'  and  the  tangent  height 
TH  is  given  in  (*).  For  a  fixed  (p,  0,  andiS, 

cos  <p'  =  cos  <p  cos  /?  -  sin  (p  sin  0  cos  0  .  (**) 

The  program  OAREJ  uses  the  values 
(p[,  RAD(l) 

<p'v  RAD<2) 

<p'v  RAD(/) 

and  maps  each  <p'  to  the  interval  [1,0]  by  means  of  the  cosine  mapping  and 
then  to  the  interval  [0,1]  using  the  map  y  =  -x.  An  interpolation  precedes 
the  last  map. 


RADINT  (J)  RADINT  tJ) 

i  M  i  i  i  i  y  x  I  i  I  i  1 1  I  i  I  l  i  i  i 


From  (**),  for  fixed  <p  and  £ 


-cos  tp'  =  -cos  <p  cos  /S  +  sin  <p  sin  0  cos  0 


(+> 


and  the  relation  (+)  is  used  as  a  map  of  0  to  -cos  <p'  and  hence  a  radiance 
value  is  associated  with  ((0,6,0). 

Measured  values  W(to.,0^)  of  the  off-axis  rejection  function  are  input 
to  OARF.J  on  cards.  For  each  fixed  tp.,  OARKJ  uses  one  of  several  options 
to  generate  interpolated  values  of  W(<p^,  *).  The  option  chosen  is  related 
to  the  number  of  0.  inputs. 

Finally,  Kq.  (G)  is  computed  as  follows: 


G(0) 


■/. 


rr/2 


-tt/2 


d0  W(tp ,  0)  I(<p 


,M)  =  zf 


it/ 2 


d0  W (<p,0)  HQ, 0,0) 


900 

-2  J^W(<p,0k)  I«p,0k,/J)  A0k 
11=0 

(the  first  and  last  terms  are  actually  modified),  where  0k  =  (k/10)  •  1°  and 
A0k  -  (1/I0)th  of  a  degree,  then 


and  so 


OAR(/3)  = 


OAR(0)  = 


rlt/ 2 

dtp  sin  tp  d9  W(tp,0)  l(<p,0,0) 


fit/ 2 

dipsintp-2l  d0  W(tp ,  0)  I(<p,  0,/3) 
J0 


ckp  sin  <p  G(<p) 


sinViG(W  +  sin^G(V 

2 


A<^ 


in  addition  to  computing  OAli(0),  the  program  modifies  OAR(/S)  using  the 
field  of  view  characteristics  of  the  sensor. 
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Subroutine  to  Extend  LOWTRAN  4  Computed  Radiance  for  Vehicle  AIH- 
tude  =  256  km  to  Other  Vehicle  Altitudes 


The  standard  input  to  the  OFF  AXIS  REJECTION  program  includes  radi¬ 
ance  values  for  25  zenith  angles  (lines  of  sight)  between  90°  and  180°  and 
for  a  vehicle  altitude  of  256  km.  If  off-axis  radiation  results  are  desired 
for  vehicle  altitudes  other  than  256  km,  then 

(a)  LOWTRAN  4  can  be  rerun  for  the  new  vehicle  altitude  and  for  appro¬ 
priate  zenith  angles. 

(b)  the  LOWTRAN  4  -  256  km  results  can  be  interpreted  as  new  vehicle 
altitude  results  by  observing  the  relationship  between  corresponding 
zenith  angles. 

An  algorithm  for  part  (b)  has  been  incorporated  in  the  OAREJ  program 
and  the  procedure  used  is  as  follows: 


The  underlying  assumption  is  that  the  radiance  incident  on  the  sensor  A 
and  on  sensor  B  along  the  indicated  line  of  sight  is  the  same  and,  in  effect, 
only  depends  on  TH. 

By  right  triangle  trigonometry. 


sin 


TH+R£ 
hTTr^  • 


sin 


TH+  Rr 
_ E_ 

H2+  Re 


which  implies 


/S'  =  arcsin 


+  Re*Re  | 

H2  *  Re  \  ' 


Consequently,  radiance  values  from  LOWTRAN  4  for  HI  =  256  km  and  angle 
/3^  will  be  the  same  as  the  radiance  for  H2  and  angle  /8  .  Thus,  if  the  fol¬ 
lowing  table  gives  the  input  to  the  off-axis  rejection  program  for  HI  =256, 


Angle 

Channel  1 

Channel  2 

Channel k 

*1 

W 

W 

•  •  • 

W 

S2 

‘A> 

'A> 

•  •  • 

"A1 

« 

• 

•  •  • 

W 

• 

K 

W 

;A» 

•  •  * 

• 

then  exactly  the  same  table  is  appropriate  for  a  vehicle  altitude  H2,  but  the 
reference  angle  /3  is  changed: 


Angle 

Channel  1 

Channel  2 

Channel  k 

W 

•  •  • 

W 

*2 

W 

•  •  • 

A> 

*3 

• 

[2(^3) 

•  •  • 

'A> 

• 

fi' 

Hn 

;a> 

w 

•  •  • 

• 

w 
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Input  to  OAREJ:  The  Off  Axis  Rejection  Curves 

For  a  given  line  of  sight  -  the  z-axis  in  the  figure  below  -  the  off-axis  re¬ 
jection  curves  are  measured  as  weights  W(<p,0).  where  o  is  the  angle  from 


the  z-axis  to  the  x-axis  and  0  is  the  angle  from  the  x-axis  to  the  y-axis. 

See  Figure  1  for  the  relation  between  the  x-,  y-,  and  z-axis  and  the 
LOWTRAN  4  coordinate  system  x',  y' ,  z' . 

Typically,  a  sequence  0  <<p0  <  •  •  •  <<pN  <;  90°  of  values  is  speci¬ 

fied.  Then,  for  each  <p.  one  or  more  0..'s,  0  <  0..  <  90°,  are  selected  and 

l  ij  ij  “ 

measurements  of  the  off-axis  weights  0..)  are  taken.  There  are  three 
options  for  the  W(cp ,  0)  input  to  OARKJ: 

1.  For  each  <p. ,  only  one  9.  is  selected  and  the  one  value  W(<p. ,  0. )  is 

i  .1  i  r 

measured.  The  assumption  is  then  made  that 

W(<0.  .9)  ,  -90°  <9  <90°  . 

2.  For  each  <0. ,  two  0?s  are  selected  0°  <  0.  <  8,  <  90rt  and  the  values 

1  .1  K 

W{(ffj,9.)  andW(<p..  0^)  arc  input.  The  assumption  in  this  case  is  that 
W«Pj,  9.)  and  ,  0^)  lie  on  an  ellipse  W{v?. ,  0).  The  parameters  de¬ 
fining  the  ellipse  are  determined  from  the  above  two  input  values. 
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3.  For  oach  <0.,  three  or  more  6.'s  are  specified.  In  this  case,  the  values 

w(<3i,e1)tw(<pi,  e2),w(<p.,  e3), . . . 


arc  used  as  a  base  to  generate  other  values  W(ex,  6)  by  interpolation. 
Note  1 :  The  assumption 

W(0,-d)  =W (<p,d)  ,  0°^dS90o 

is  in  force  throughout  OAREJ. 

Note  2:  For  Option  1  above,  the  measured  value 

w«p.,e.) 

is  actually  input  twice  and  then  Option  2  is  executed. 
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POWER  SPECTRUM  ANALYSIS,  DIGITAL  FILTERING 


Initiator:  E.A.  Murphy 

Problem  No. :  4954 
Project  No. :  2310 

This  problem  concerned  experimental  data  with  extremely  poor  signal -to- 
noise  ratio.  The  objective  was  to  investigate  the  use  of  digital  signal 
processing  techniques  to  improve  the  signal/noise  characteristic  in  order 
to  extract  desired  information. 

A  program  was  developed  to  examine  the  frequency  spectra  of  the 
data. 

This  program  calculates  the  Fast  Fourier  Transform  of  the  specified 
input  curve.  It  prints  its  magnitude  and  phase  and  then  plots  its  normal¬ 
ized  magnitude.  The  core  storage  needed  (as  specified  in  program)  is 
directly  related  to  the  number  of  Fast  Fourier  points  desired.  The  less 
points  desired  the  less  core  storage  required. 

Coupled  with  the  high  noise  level,  it  was  found  that  the  data  drop  out 
rate  was  so  high  that  severe  aliasing  of  the  spectra  was  present.  This 
rendered  the  extractions  of  useful  information  hopeless.  Thus,  no  posi¬ 
tive  results  could  be  obtained . 
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AUTOMATED  AZIMUTH  MONITORING 


Initiator:  T.E.  Wirtanen 

Problem  No.:  4997 
Project  No. :  7600 

The  purpose  of  this  problem  was  to  extract  data  from  a  computer  tape 
which  was  generated  by  an  "in-house"  automated  azimuth  monitoring 
system. 

After  runni ng  several  tape  dumping  programs,  it  became  apparent 
that  most  of  what  was  occurring  was  meaningless  noise,  because  all  the 
meaningful  information  was  packed  into  the  very  first  record,  which  was 
not  getting  picked  up  correctly. 

As  it  turned  out,  the  mistakes  are  occurring  at  a  hardware  level, 
probably  because  the  designers  of  the  tape  controller  were  unaware  of 
industry-wide  standards  used  by  tape  transport  manufacturers.  This  is 
much  more  fundamental  than,  for  example,  reading  a  tape  written  in  a 
"wrong"  or  "foreign"  format. 

Thus,  the  problem  was  returned  to  the  initiator  with  some  information 
as  how  to  approach  it. 


112 


